We  present  a  new  nonlinear  optimization  procedure  for  the  computation  of  generalized 
Gaussian  quadratures  for  a  broad  class  of  functions.  While  some  of  the  components 
of  this  algorithm  have  been  previously  published,  we  present  a  simple  and  robust 
scheme  for  the  determination  of  a  sparse  solution  to  an  underdetermined  nonlinear 
optimization  problem  which  replaces  the  continuation  scheme  of  the  previously  pub¬ 
lished  works.  The  performance  of  the  resulting  procedure  is  illustrated  with  several 
numerical  examples. 
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including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

30  JUN  2008  2- REPORT  TYPE 

3.  DATES  COVERED 

00-00-2008  to  00-00-2008 

4.  TITLE  AND  SUBTITLE 

A  Nonlinear  Optimization  Procedure  for  Generalized  Gaussian 
Quadratures 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROTECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Yale  University  , Department  of  Computer  Science, New 

Haven, CT, 06520-8285 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBIECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

_ _ _  ABSTRACT 

18.  NUMBER  19a.  NAME  OF 

OF  PAGES  RESPONSIBLE  PERSON 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  Same  OS 

unclassified  unclassified  unclassified  Report  (SAR) 

34 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


2 


1.  Introduction 


Classical  Gaussian  quadrature  rules  are  optimal  formulas  for  the  evaluation  of 
integrals  of  the  form 

(1.1)  f  P(x)u(x)dx 

J  a 

where  P(x)  is  a  polynomial  and  u  :  [a,  b\  — »  M+  is  an  essentially  arbitrary  nonnegative 
weight  function;  they  are  optimal  in  sense  that  an  n-point  Gaussian  quadrature  rule 
integrates  polynomials  of  degree  n  —  1  exactly  with  respect  to  the  weight  function 
u.  They  are  extremely  efficient  as  long  as  functions  to  be  integrated  are  smooth  or 
of  the  form  f(x)  =  g(x)w(x)  where  g(x)  is  smooth  and  w(x)  is  fixed  and  known  a 
priori.  However,  they  perform  poorly  in  a  number  of  instances  of  great  importance 
in  computational  mathematics.  Of  particular  interest  is  the  integration  of  functions 
/  which  admit  representations  of  the  form 


(1.2)  f(x)  =  ^2ai(j>i(x), 

1=1 


where  the  0j  are  oscillatory,  singular,  or  both  and  known  a  priori  but  the  coefficients 
at  are  not.  Then  each  of  the  functions  (f>i(x)  can  exhibit  a  different  type  of  singularity 
or  a  different  frequency  of  oscillation,  and  classical  Gaussian  quadrature  rules  perform 
very  poorly.  Many  such  integrals  arise  in  computational  physics  (see,  for  instance, 
[4]  and  [18]). 

It  has  long  been  known  that  Gaussian  quadratures  admit  generalization  to  a  fairly 
broad  class  of  systems  of  functions  (see  [11,  12,  15,  16,  13]).  We  adopt  the  termi¬ 
nology  of  [14,  20,  3],  and  say  that  a  generalized  Gaussian  quadrature  for  a  collection 
{0i,...,02  n}  of  2  n  functions  and  a  nonnegative  weight  function  oj  is  an  n-point 
quadrature  rule 


(1.3) 


n 

E 

3= 1 


< Kxi 


)Wi 


4>{x)uj{x)dx 


exactly  integrating  each  of  the  0,  with  respect  to  the  weight  function  u. 

The  constructions  found  in  [11,  12,  15,  16,  13]  do  not  lead  readily  to  numerical  al¬ 
gorithms  for  the  design  of  generalized  Gaussian  quadrature  rules.  In  [14] ,  a  numerical 
algorithm  is  introduced  for  the  construction  of  generalized  Gaussian  quadrature  rules 
for  a  fairly  broad  class  of  functions.  The  approach  is  based  on  the  observation  that 
the  nodes  x3  and  weights  Wj  of  a  generalized  Gaussian  quadrature  satisfy  a  nonlinear 
system  of  equations.  The  procedure  of  [14]  is  a  variant  of  Newton’s  algorithm  cou¬ 
pled  with  a  continuation  scheme  for  the  generation  of  a  suitable  initial  point  for  the 
modified  Newton  iterations.  In  [20]  and  [3],  the  continuation  scheme  of  [14]  is  refined, 
improving  the  stability  of  the  algorithm,  and  a  new  preprocessing  step  is  added  in 
[20],  greatly  expanding  its  range  of  applicability. 

The  present  paper  introduces  a  new  numerical  procedure  for  obtaining  generalized 
Gaussian  quadrature  rules.  While  we  also  use  Newton-type  iterations  to  solve  nonlin¬ 
ear  systems  of  equations,  our  scheme  differs  from  the  algorithm  of  [14,  20,  3]  in  that 
the  continuation  method  is  abandoned  in  favor  of  a  simpler,  more  robust  approach. 


3 


The  paper  is  structured  as  follows.  Section  2  contains  mathematical  and  numerical 
preliminaries.  In  section  3,  we  describe  certain  standard  numerical  tools  used  by  the 
algorithm.  Section  4  describes  the  algorithm  in  detail.  Section  5  contains  several  ex¬ 
amples  of  quadratures  we  have  obtained.  Finally,  in  Section  6,  we  present  conclusions 
and  discuss  several  possible  extensions  of  this  work. 


2.  Mathematical  and  numerical  preliminaries 

2.1.  Generalized  Gaussian  and  Chebyshev  quadratures.  The  quadrature  for¬ 
mulas  considered  in  this  paper  are  of  the  form 

n 

(2.1) 

3= 1 

where  the  Xj  and  c Oj  are  real  numbers.  We  will  refer  to  the  Xj  as  the  nodes  and  the 
Wj  as  the  weights  of  the  quadrature  formula  (2.1).  They  will  be  used  to  approximate 
integrals  of  the  form 

(2.2)  f  (p(x)u(x)dx 

J  a 

where  uj{x)  is  a  nonnegative  weight  function. 

Quadratures  are  typically  chosen  so  as  to  make  (2.1)  exact  for  some  set  of  functions, 
commonly  polynomials  of  a  fixed  order.  Classical  Gaussian  quadrature  rules  consist 
of  n  nodes  and  n  weights  which  integrate  polynomials  of  order  2n  —  1  exactly.  In  [14], 
the  notion  of  a  Gaussian  quadrature  was  generalized  as  follows: 

DEFINITION  2.1.  A  quadrature  formula  will  be  referred  to  as  Gaussian  with  re¬ 
spect  to  a  set  of ‘In  functions  (pi, ... ,  (p2n  :  [a,  b]  — ■>  M  and  a  weight  function  oj  :  [a,  b]  — > 
M+  if  it  consists  of  n  nodes  and  n  weights  and  integrates  exactly  the  functions  (pi  with 
respect  to  the  weight  function  u  for  all  i  =  1, . . . ,  2 n.  The  weights  and  nodes  of  a 
Gaussian  quadrature  will  be  referred  to  as  Gaussian  weights  and  nodes,  respectively. 

Although  Chebyshev  quadratures  are  classical  Gaussian  quadratures  on  the  interval 
[—1, 1]  with  respect  to  the  weight  function  u(x)  =  (1  —  x2)~1^2,  in  practice,  Chebyshev 
nodes  and  weights  are  often  used  to  integrate  functions  on  [—1,1]  with  respect  to 
the  weight  function  uj(x)  =  1.  This  practice  leads  to  a  2n-point  quadrature  which 
integrates  exactly  polynomials  of  order  2n  —  1,  and  motivates  the  following  definition: 

DEFINITION  2.2.  A  quadrature  formula  will  be  referred  to  as  a  Generalized  Cheby¬ 
shev  quadrature  for  a  set  of  2n  functions  (pi, ... ,  (p2n  ■  [a,  b]  — >  R  and  a  weight  function 
oj  :  [a,  b]  — >  M+  if  it  consists  of  2 n  nodes  and  2 n  weights  and  integrates  exactly  the 
functions  (pi  with  respect  to  the  weight  function  oj  for  all  i  =  1, . . ,  2 n.  The  weights 
and  nodes  of  a  Chebyshev  quadrature  will  be  referred  to  as  Chebyshev  weights  and 
nodes,  respectively. 


2.2.  Quadrature  rules  and  optimization.  Let  m,n  be  integers  with  m  <  2 n  and 

let  aq,,* . . ,  xn,  W\, . . . ,  wn  be  a  quadrature  rule 


x)u{x)dx  ~  ^  f(xi) 

3= 1 


un 


(2.3) 
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which  is  exact  for  the  functions  0i Obviously,  the  weights  Wj  and  nodes 
Xj  of  such  a  quadrature  satisfy  the  (generally  underdetermined)  nonlinear  system  of 


equations 

Fi(xi,  .  . 

•  i  )  •  •  • 

,  Wn)  =  bi 

(2.4) 

F2(x i, . . 

•  >  7  •  •  • 

,  Wn)  =  b2 

Fm  (*£l  i  ■  ■ 

•  7  7  •  •  • 

i  ®n)  bm , 

where 

(2.5) 

Fi(x  i,  ...,xn, 

»'i,  ■••,»«) 

n 

=  ^HX3 

3= 1 

and 

(2.6) 

bi  = 

:  /  <f>i(x)u( 

J  a 

[x)dx. 

When  m  =  2 n,  the  nodes  x3  and  weights  Wj  in  (2.4)  form  a  generalized  Gauss¬ 
ian  quadrature  for  the  functions  <f>i, . . . ,  02n-  In  [11]  and  [13],  the  existence  of  a 
unique  solution  of  (2.4)  is  proven  under  certain  conditions  on  the  system  of  functions 
0i, ... ,  02 n-  It  is  observed  in  [14,  20,  3]  that  solutions,  or  at  least  approximate  solu¬ 
tions,  of  (2.4)  exist  for  many  systems  of  2 n  functions  not  satisfying  those  conditions. 

In  this  paper,  we  will  encounter  systems  of  the  form  (2.4)  where  m  <  2 n.  Under 
these  conditions,  the  system  does  not  admit  a  unique  solution.  Non-uniqueness  is  not, 
however,  an  obstruction  to  the  determination  of  a  quadrature  rule  for  the  functions 
0i, ... ,  0m  since  such  a  rule  is  given  by  any  solution  of  the  system  (2.4).  Indeed,  even 
in  the  case  when  there  is  no  exact  solution  of  (2.4),  it  is  often  possible  to  construct 
an  approximate  quadrature  for  the  functions  0i, . . . ,  0m. 

We  close  this  section  with  the  following  definition: 

DEFINITION  2.3.  Suppose  0i,...,0m  :  [a,  b]  — >  M  are  square  integrable  with  re¬ 
spect  to  the  nonnegative  integrable  weight  function  u.  We  say  that  a  quadrature  rule 


Xi,  •  •  •  ,£n,Wl,  . 

. . ,  vjn  integrates  0i, . . . ,  0m 

with  precision  e  >  0  if 

(2.7) 

m 

\Fi(xi, . . xn.  «•,. 

2=1 

■  ■  ■  ,wn)  -  bf2  <  e2, 

where 

(2.8) 

Fi(x i,  •  •  •  ,xn,wu  . . . . 

n 

,Wn)  =  ^2<j>i(Xj)Wj 

3= 1 

and 

(2.9) 

bi  =  00 

J  a 

x)u>(x)dx 

for  all 

. ,  m. 
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2.3.  Quadrature  and  interpolation.  It  is  well  known  that  when  Chebyschev  nodes 
are  used  for  polynomial  interpolation  on  the  interval  [—1,1],  the  procedure  is  numer¬ 
ically  stable  and  the  convergence  properties  are  close  to  optimal  (see  [19]  and  [6]). 
In  this  subsection,  we  prove  that  the  nodes  of  any  Gaussian  quadrature  and  many 
generalized  Gaussian  quadratures  lead  to  stable  interpolation  formulas. 

The  principal  analytic  tool  of  this  subsection  is  the  following  obvious  theorem  (see, 
for  example,  [20]). 

THEOREM  2.1.  Suppose  that  the  weight  function  u>  :  [a,  b]  — >  M  is  nonnegative 
and  the  functions  (pi, ...  ,(pn  :  [a,b]  — ■>  R  are  orthonormal  with  respect  to  u.  Suppose 
further  that  the  n-point  quadrature  rule  x\, ... ,  xn,  w±, . . . ,  wn  is  exact  for  all  products 
of  the  form  (pi(x)cpj(x)  and  Wi  >  0  for  all  1  <  i  <  n.  Then  the  n  x  n  matrix  A  with 
entries 

(2.10)  Ai:j  =  ^fwf(pi(xj) 
is  orthogonal. 

Now  let  /  be  a  function  of  the  form 

n 

(2.11)  f{x)  =  ^2,aif)i{x). 

3= 1 

We  would  like  to  construct  an  interpolation  formula  on  the  interval  [a,b]  for  functions 
of  this  form;  that  is,  given  the  values  fi, . . . ,  fn  of  a  function  of  the  form  (2.11) 
at  a  collection  of  points  x\,...,xn,  we  would  like  a  formula  for  determining  the 
coefficients  an, . . . ,  an.  Let  a  be  the  vector  a  =  (an, . . . ,  an),  let  F  be  the  vector 
F  =  (/i,  /2, . . . ,  fn),  and  finally,  let  B  be  the  n  x  n  matrix  with  entries 

(2.12)  Bij  =  faixj). 

Then 

(2.13)  F  =  Ba, 

and  assuming  that  B  is  invertible,  it  follows  that 

(2.14)  a  =  B^F. 

If  the  condition  number  of  B  is  not  too  high,  then  (2.14)  is  a  numerically  stable 
formula  for  computing  the  coefficients  an,.-.. ,  an.  The  following  observation  is  the 
principal  observation  of  this  subsection. 

Observation  2.1.  Under  the  conditions  of  Theorem  2.1,  the  matrix  B  in  the  inter¬ 
polation  formula  (2.1 4)  is  of  the  form 

(2.15)  B  =  DQ, 

where  D  is  a  diagonal  matrix  with  entries 

(2.16)  Dn  =  -L 

Vwi 

and  Q  is  orthogonal.  Thus 

(2.17)  a  =  Q*D~1F. 

In  other  words,  the  coefficients  a  can  be  obtained  by  applying  a  diagonal  matrix 
followed  by  an  orthogonal  matrix. 
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2.4.  The  damped  Gauss-Newton  method.  The  damped  Gauss-Newton  method 
is  a  well-known  iterative  technique  for  the  solution  of  nonlinear  least-squares  prob¬ 
lems.  It  converges  under  very  general  conditions,  and  does  not  require  that  the 
Jacobian  of  the  system  be  nonsingular.  Here  we  give  only  elementary  details;  a  more 
thorough  treatment  can  be  found  in  [7]. 

Suppose  that  R  :  Rn  — »•  Mm  is  a  continuously  differentiable  function  of  the  form 


(2.18) 


R(x) 


/  n(x)  \ 

r2  (x) 


V  rm(x)  ) 


and  let  J(x)  be  the  Jacobian 


(2.19) 


J(x) 


few 

few 


few  \ 
few  / 


of  R  at  the  point  x.  The  damped  Gauss-Newton  method  is  a  numerical  method  for 
minimizing  the  function 


(2.20) 


1  m 

f(x)  =  2  \rj(x)\2- 

3= 1 


It  belongs  to  a  broad  class  of  numerical  optimization  methods  which  proceed  from 
an  initial  guess  x0  by  forming  a  sequence  x±,  x2,  ■  ■  ■  of  iterates  via  the  formula 


(2.21) 


xk+ 1  xk  T 


where  dk  is  referred  to  as  the  search  direction  at  iteration  k  and  a k  is  a  carefully 
chosen  step  size. 

The  primary  purpose  of  this  section  is  Theorem  2.2  below,  which  gives  conditions 
under  which  an  iteration  of  the  form  (2.21)  converges.  We  start  with  the  following 
definition: 

DEFINITION  2.4.  Let  f  :  Wl  M  be  a  continuously  differentiable  function  and 
consider  any  iteration  of  the  form  (2.21).  We  say  that  the  step  length  ak  and  step 
direction  dk  satisfy  the  Wolfe  conditions  at  the  point  Xk  if 


(2.22)  f(xk  +  akdk )  <  f(xk)  +  \akW f(xk)  ■  dk, 
and 

(2.23)  V/(xfc  +  akdk)  ■  dk  >  fN f{xk)  ■  dk 


for  some  constants  A  G  (0, 1)  and  (3  G  (A,  1). 

The  following  theorem  can  be  found  in  [7]: 

THEOREM  2.2.  Suppose  that  f  :  Mn  — >  M  is  a  continuously  differentiable  function 
bounded  from  below,  and  assume  that  there  exists  7  >  0  such  that 

(2.24) 


l|V/(x)  —  V/(y)||2  <  ^\\x  —  y\\2 
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for  every  x  and  y  inW1 .  Consider  any  iteration  of  the  form  (2.21)  such  that  for  each 
k  =  0, 1, . . .  the  pair  (dk,  o>k)  satisfies  the  Wolfe  conditions.  If,  in  addition,  one  of  the 


conditions 

(2.25) 

V/(xfc)  •  dk  <  0 

or 

(2.26) 

V f(xk)  =  0  and  dk  =  0 

holds  for  each  k  =  0, 1, . . 

.,  then  either 

(2.27) 

V f(xk)  =  0  for  some  k  >  0, 

or 

(2.28) 

,.  Vf{xk)-dk 
hm  ,,  ,  ,,  =  0. 

fc^oo  1 1  Gtfc  1 1 2 

Remark  2.1.  Theorem  2.2  states  that  either  the  sequence  Xk  converges  to  a  critical 
point  for  the  function  f  or  the  direction  dk  become  orthogonal  to  the  gradient  of  f . 
In  practice,  it  is  easy  to  avoid  the  later  condition,  thus  ensuring  the  convergence  of 
V/(xfc)  to  0. 

Remark  2.2.  Under  mild  conditions  on  the  function  f ,  given  a  sequence  of  {dk}  sat¬ 
isfying  condition  (2.25),  there  exists  a  sequence  of  ak  satisfying  the  Wolfe  conditions 
(see,  for  example,  [7]  Theorem  6.3.2). 

In  the  case  of  the  damped  Gauss-Newton  method,  the  search  direction  dk  is  chosen 
as  a  solution  to  the  least-squares  problem 

(2.29)  min  |[ J(xk)dk  +  R(xk)\\2, 

dk 

which  is  an  affine  approximation  to  the  nonlinear  least-squares  problem 

min  ||-R(x)||2. 


Since 

(2.30) 

V/(xfc)  =  J*(xk)R(xk), 

we  have 

(2.31) 

(V/(xfc),4)  =  (R(xk),J(xk)dk) 

The  search  direction  dk  is  chosen  so  that  J(xk)dk  is  the  projection  of  —R(xk)  onto 
the  column  space  of  J(xk).  It  follows  that 

(2.32)  (Vf(xk),dk)  =  (R(xk),  J(xk)dk)  <  0. 

If  we  choose  dk  to  be  0  in  the  event  that  (V f(xk),  dk)  =  0,  then  we  obtain  the 
following  theorem: 

THEOREM  2.3.  Suppose  that  R  :  Rn  — >  Mm  is  a  continuously  differentiable  func¬ 
tion  of  the  form  (2.18),  f  :  Mn  — »•  R  is  given  by  (2.20),  and  the  Jacobian,  J(x),  of  R 
is  given  by  (2.19).  Further  suppose  that  ||  J(x)\\2  is  bounded  on  Rn  and  there  exists  a 
constant  7  >  0  such  that 


\\J{x)-J{y)h<7\\x-y\\2 
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for  all  x  and  y  in  Mn.  If  xk  is  defined  by  the  iteration 

Xk+ 1  “1“  Olkdk 

where,  for  each  k  =  1,  2, . . dk  is  a  solution  of  the  least-squares  problem 
II J{xk)dk  +  R{xk) II2  =  min  \\J{xk)v  +  R(xk) ||2, 

V 

and  the  sequence  {ak}  is  chosen  so  that  the  pairs  (dk,  oik)  satisfy  the  Wolfe  conditions 
for  k  =  0, 1, . . then  either 

V f(xk)  =  0  for  some  k  >  0 


or 


lim 

k — »oo 


YfM ' dk 

\\dk  ||  2 


0. 


2.5.  Singular  value  decomposition.  The  SVD  is  a  ubiquitous  tool  in  numerical 

analysis  (see,  for  instance,  [8]).  Here  we  discuss  it  in  the  case  of  real  matrices. 

LEMMA  2.1.  (SVD).  For  any  n  x  m  real  matrix  A,  there  exist,  for  some  integer 
k,  an  n  x  k  real  matrix  U  with  orthonormal  columns,  an  m  x  k  real  matrix  V  with 
orthonormal  columns,  and  akxk  real  diagonal  matrix  £  with  positive  diagonal  entries 
Gy  >  a 2  >  •  •  •  >  (Tk  >  0,  such  that 


(2.33) 


A  =  U  -T,  -V*. 


The  diagonal  entries  of  £  are  called  singular  values,  the  columns  of  the  matrix  U 
are  called  the  left  singular  vectors,  and  the  columns  of  the  matrix  V  are  called  right 
singular  vectors. 

One  of  the  most  common  applications  of  the  SVD  is  the  approximation  of  matrices; 
if  we  let  £p  denote  the  diagonal  matrix  with  entries 


then 

(2.34) 


Dr 


Oi  if  i  <  p 
0  otherwise  ’ 


\\A-UFpV*\\2  =  ap+l. 


The  matrix  UY,pV*  is,  in  fact,  the  optimal  rank  p  approximation  of  the  matrix  A 
(see,  for  instance,  [8]);  that  is, 

(2.35)  min  \\A  -  Ak ||2  =  ap+l 


where  Ak  ranges  over  the  set  of  all  n  x  m  matrices  of  rank  k. 


2.6.  QR  decompositions.  The  singular  value  decomposition  provides  the  means  to 
construct  an  optimal  rank  k  approximation  to  a  given  matrix;  however,  the  SVD  can 
be  expensive  to  form,  and  other  less  computationally  expensive  matrix  factorizations 
can  be  used  to  compress  matrices  in  lieu  of  the  SVD. 

By  a  slight  abuse  of  terminology,  we  will  refer  to  the  factorization  in  the  following 
obvious  lemma  as  a  “QR  decomposition.” 
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LEMMA  2.2.  For  any  n  x  m  real  matrix  A,  there  exist  an  integer  k,  an  n  x  k  real 
matrix  Q  with  orthonormal  columns,  an  m  x  m  permutation  matrix  II,  and  a  k  x  m 
real  matrix  R  with  the  block  form 

R=  (T  M  )  , 

where  T  is  a  k  x  k  upper  triangular  matrix  with  positive  diagonal  entries  and  M  is  a 
general  k  x  (m  —  k)  matrix,  such  that 

(2.36)  An  =  QR. 

Moreover,  there  is  a  one-to-one  correspondence  between  m  x  m  permutation  matrices 
n  and  decompositions  of  the  form  (2.36). 

The  following  theorem  can  be  found  (in  a  stronger  form)  in  [10] . 

THEOREM  2.4.  Suppose  that  A  is  a  real  m  x  n  matrix  with  singular  values  oq  > 
O2  >  •  •  •  >  ar  >  0.  For  any  integer  p,  0  <  p  <  r,  there  exists  an  m  x  m  permutation 
matrix  n  such  that  the  QR  decomposition  uniquely  determined  by  n  is  of  the  form 

(2.37) 

where  Rn  is  a  p  x  p  upper  triangular  matrix  with  positive  diagonal  entries,  and 

(2.38)  H-R22II2  <  \A  +  (n  -  k)crp+i. 


Remark  2.3.  Theorem  2. 4  implies  that  given  any  real  m  x  n  matrix  A  with  singular 
values  oq  >  cr2  >  . . .  >  <rr  >  0  and  any  integer  0  <  k  <  r,  there  is  a  permutation  ma¬ 
trix  n  such  that  the  well-known  Gram- Schmidt  orthogonalization  procedure  (see,  for 
example,  [8]j  applied  to  the  matrix  All  results  in  the  approximate  matrix  factorization 

(2.39)  \\A-QS\\2<cak+1, 

where  Q  is  an  mxk  matrix  with  orthonormal  columns  and  S  is  a  kxn  matrix  S.  Note 
that  the  error  bound  (2.39)  is  close  to  the  optimal  error  for  rank  k  approximations  to 
A. 

Remark  2.4.  When  properly  implemented,  the  modified  Gram-Schmidt  algorithm 
with  reorthogonalization  is  a  reliable  method  for  obtaining  a  QR  decomposition  (see, 
for  instance,  [2]).  While  there  are  no  guaranteed  bounds  of  the  form  (2.39)  for  this 
algorithm,  it  does  well  in  practice.  In  [10],  robust,  provably  stable  algorithms  are 
introduced  for  producing  QR  decompositions  that  satisfy  bounds  of  the  form  (2.39). 

2.7.  Analogs  of  matrix  factorizations  for  finite  sequences  of  functions.  Here 
we  introduce  analogs  of  the  matrix  factorizations  of  the  preceding  two  subsections 
for  finite  sequences  of  functions.  We  begin  with  the  following  obvious  generalization 
of  Lemma  2.2. 

LEMMA  2.3.  For  any  finite  sequence  <j)  1, . . . ,  (pm  :  [a,  b]  — >  M  of  square  integrable 
functions,  there  exists  an  integer  k  <  m,  a  permutation  n  e  Sm,  a  collection  of 
orthonormal  functions  qi, . . . ,  :  [a,  b]  — >  M,  and  an  k  x  m  real  matrix  R  =  (rtj) 
with  the  block  form 


(  T  M  ) 
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where  T  is  a  k  x  k  upper  triangular  matrix  with  positive  diagonal  entries  mid  M  is  a 
general  k  x  (rn  —  k )  matrix,  such  that 

max(jf,/c) 

(2.40)  </>n  (j)(x)  = 

i=  1 

for  all  j  =  1 ,m.  Moreover,  there  is  a  one-to-one  correspondence  between  decom¬ 
positions  of  this  form  and  permutations  II  G  Sm. 

By  analogy  with  the  finite  dimensional  case,  we  will  refer  to  decompositions  of 
this  type  as  QR  decompositions.  The  functions  (pix)  will  be  referred  to  as  QR 
functions  and  we  will  call  the  diagonal  entries  rlt  of  R  the  normalizing  factors  of  the 
decomposition. 

As  in  the  case  of  matrices,  a  common  application  of  expansions  of  the  form  (2.40) 
is  the  “compression”  of  the  functions  (bp  that  is,  often  permutations  II  G  Sm  can  be 
found  for  which  the  truncated  series 

p 

(2.41)  ^OijQiix), 

3= 1 

with  p  <  k,  provide  good  approximations  to  the  functions  4>n(j)- 

Now,  we  restate  a  result  found  in  [3],  generalizing  the  SVD  to  the  case  of  a  finite 
sequence  of  functions. 

THEOREM  2.5.  Suppose  that  the  functions  (f>i,  '■  [«,  b]  — >  M  are  square 

integrable.  Then  there  exist  an  integer  k,  a  finite  orthonormal  sequence  of  functions 
Mi, . . . ,  Uk  :  [a,  b]  — >  M,  an  m  x  k  matrix  V  =  (ty,)  with  orthonormal  columns,  and  a 
sequence  Si  >  s2  >  •  •  •  >  Sk  >  0  G  M  such  that 

k 

(2.42)  ( bj{x )  =  y ^Ui^SjVji 

i—  1 

for  all  x  G  [a,  b]  and  all  j  =  1 , ,m.  The  sequence  s±, ...  ,Sk  is  uniquely  determined 
by  k. 

By  analogy  with  the  finite-dimensional  case,  we  will  refer  to  this  decomposition 
as  the  SVD  of  a  finite  sequence  of  functions.  We  call  the  functions  Ui  the  singular 
functions,  the  columns  of  V  the  singular  vectors,  and  the  values  s*  the  singular  values. 
The  SVD  is  clearly  a  useful  tool  for  the  compression  of  the  sequence  (pi, ,  0m:  if 
we  let  4>j{x)  denote  the  p-terrn  truncation 

p 

(2.43)  j(x )  =  y'ui(x)sivji 

i=  1 

of  the  sum  (2.42),  then 

(2.44)  \\4>j(x)  -  (bj(x)\\  <  sp+1 
for  j  =  1, . . . ,  m. 

However  it  is  obtained,  using  an  approximate  representation 

p 

f>j(x)  ~  ^2<Xiqi(x), 
i= 1 


(2.45) 
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integrals  of  the  form 
(2.46) 

can  be  approximated  as 


4>j{x)uj{x)dx 


-b  /  p 


4>j{x)oj(x)dx  ~  £  ctiqi(x)  I  uj{x)dx 


„  1=1 


(2.47) 


=  X/  a%  I  qi(x)u(x)dx-, 


i=  1 


thus,  a  quadrature  which  is  exact  for  the  integrals 


(2.48) 


qi(x)uj(x)dx 


for  i  —  1 , . . . .  ]>.  can  be  used  as  an  approximate  quadrature  for  the  integrals 


(2.49) 

for  j  —  1, . . . ,  m. 


4>j{x)ui{x)dx 


2.8.  Underdetermined  systems.  The  purpose  of  this  short  subsection  is  the  state¬ 
ment  of  two  lemmas  on  the  existence  and  behavior  of  sparse  solutions  of  underdeter¬ 
mined  systems. 

We  begin  with  the  following  result  on  the  existence  of  well-behaved  solutions  to 
underdetermined  least-squares  problems,  whose  proof  is  an  elementary  exercise  in 
linear  algebra. 

LEMMA  2.4.  Suppose  that  A  is  an  n  x  m  matrix,  n  <  m,  with  left  singular  vectors 
ui, . . .  ,Uk  and  singular  values  o\  >  . . .  >  cq,  >  0.  For  0  <  p  <  k,  let  Up  be  the 
subspace  ofWn  spanned  by  u\,. .. ,  up,  and  let  Projp  :  — >  Up  denote  the  projection 

operator  Mn  — >  Up.  Then,  given  any  vector  b  e  and  any  integer  0  <  p  <  k,  there 
exists  a  vector  x  G  with  no  more  than  p  nonzero  entries  such  that 

(2.50)  \\Ax  -  b\\2  =  \\Projpb  -  b\\2 
and 

(2.51)  \\x\\2  < —\\b\\2. 

(Jp 


Let  A  be  an  n  x  m  matrix  with  n  <  m  and  consider  the  least-squares  problem 
(2.52)  min  \\Ax  ^  b\\2, 

X 

where  b  is  a  given  vector  in  Rn.  The  problem  is  underdetermined  and  therefore  the 
condition  number  of  A  is  infinite.  However,  Lemma  2.4  implies  that  if  b  is  in  the  span 
of  a  small  number  of  singular  vectors  of  A,  then  a  well-behaved  approximate  solution 
x  to  the  least-squares  problem  (2.52)  can  still  be  found. 

The  second  lemma  is  a  stronger,  but  more  specialized  result,  whose  proof  is  an 
easy  corollary  of  the  following  theorem,  which  appears  as  Theorem  2  in  [17]. 
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THEOREM  2.6.  Suppose  that  S  is  an  arbitrary  set,  n  is  a  positive  integer,  fi, . . . ,  fn 
are  bounded  complex-valued  functions  on  S,  and  e  is  a  positive  real  number  such  that 

(2.53)  e  <  1. 

Then,  there  exist  n  points  xi, ...  ,xn  in  S  and  n  functions  gi, ...  ,gn  on  S  such  that 

(2.54)  \gk(x)\  <1  +  6 
for  all  x  in  S  and  k  —  1,  2, . . . ,  n,  and 

n 

(2-55)  f(x)  =  ^2f(xk)gk(x) 

k= 1 

for  all  x  in  S  and  any  function  f  defined  on  S  via  the  formula 

n 

(2.56)  f{x)  =  ^2,ckfk(x). 

k=  1 

Moreover,  if  the  set  S  is  finite,  then  gi, ...  ,gn  can  be  chosen  so  that  (2.5 f)  holds  with 
e  =  0. 

LEMMA  2.5.  If 

(2.57)  Ax  =  b, 

where  A  is  an  m  x  n  matrix  of  rank  m,  then  there  exists  a  vector  x  G  with  no 
more  then  m  nonzero  entries  such  that 

Ax  =  b , 

and 

||5||i  <  m||x||i. 


Proof:  By  Theorem  2.6,  there  exists  an  m  x  n  matrix  G  whose  entries  are  bounded 
by  1  and  an  m  x  m  matrix  of  A  consisting  of  m  columns  An , ... ,  Aim  of  A  such  that 

(2.58)  A  =  AG. 

Since  Ax  =  b ,  it  follows  that 

(2.59)  AGx  =  b. 

If  we  let  y  =  Gx  G  Mm  and  define  x  by  the  formula 


Xj  = 


(2.60) 


Vk  if  j  =  h, 

0  otherwise  ' 
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then  x  has  at  most  m  nonzero  entries,  Ax  =  b ,  and 

m 

Piii  =  ^2  1^1 


(2.61) 

(2.62) 

(2.63) 

(2.64) 


i— 1 


m 


E 

i— 1 


n 


9ijXj 


3  = 1 


*=i  j=i 
<  m||x||i, 


i, 


as  desired. 

Remark  2.5.  For  sets  S  which  are  finite,  the  proof  of  Theorem  2.6  given  in  [17] 
is  constructive,  but  the  procedure  is  computationally  infeasible.  To  wit,  in  the  spe¬ 
cial  case  of  Lemma  2.5,  the  scheme  of  [17]  amounts  to  the  choice  of  a  submatrix  A 
consisting  of  a  set  of  columns  Aix , ... ,  Aim  of  A  which  maximize  the  determinant 

(2.65)  det  (  Ah  Ah  ...  Aim  ) 

over  the  collection  of  all  sets  of  m  columns  of  the  matrix  A.  The  existence  of  the 
matrix  G  follows  since 

det  (A)  0 

by  construction,  and  the  bound  on  the  entries  of  G  follows  from  Cramer’s  rule. 

Remark  2.6.  In  practice,  the  modified  Gram-Schmidt  procedure  with  reorthogonal- 
ization  can  be  used  to  find  a  set  of  m  columns  An , ... ,  Aim  of  the  m  x  n  rank  m 
matrix  A  such  that 


(2.66)  det  (  Ah  Ai2  ...  Aim  ) 
is  comparable  to  the  supremum 

(2.67)  sup  det  (  An  Ah  ...  Ajm  ) 

over  all  collections  of  m  columns  of  A.  Indeed,  the  modified  Gram-Schmidt  procedure 
is  nothing  more  than  a  “greedy”  algorithm  for  finding  an  approximate  solution  to  the 
optimization  problem 

(2.68)  argmax  det  (  An  Aj2  . . .  Ajm  )  . 

ilr-jm 

An  obvious  modification  of  the  argument  in  [17]  (which  is  sketched  above  in  Remark 
2.5)  shows  that  once  such  a  set  of  columns  i\, . . .  ,im  has  been  found,  the  matrix  A 
can  factored  as 

A  =  AG, 

where  A  is  an  m  x  m  submatrix  of  A  and  G  is  an  m  x  n  matrix  whose  entries  are 
bounded  in  absolute  value,  but  not  necessarily  by  1. 

Remark  2.7.  Strong  rank  revealing  QR  factorizations  (see  [10]  and  [2]  ,  for  instance) 
identify  a  spanning  set  An , ... ,  Aim  of  columns  of  an  m  x  n  matrix  A  of  rank  m  for 
which  the  ratio  of 


det  (  An  Ai2  ...  Aim  ) 


14 


to 


sup  det  (  Ah  Aj2  . . .  Ajm  ) 

ilv-jm 

is  guaranteed  to  satisfy  a  lower  bound,  thus  ensuring  that  a  stable  factorization  of  the 
form 

A  =  AG, 


where  A  is  an  m  x  m  submatrix  of  A,  can  be  found. 


2.9.  The  Sherman-Morrison- Woodbury  formula.  The  Sherman-Morrison- Woodbury 
formula  gives  an  expression  for  the  rank-k  update 

(A  +  UVY1 

of  the  inverse  of  a  matrix  A.  The  following  Lemma  can  be  found,  for  instance,  in  [8]: 

LEMMA  2.6.  Suppose  that  A  is  an  invertible  n  x  n  matrix,  U  is  an  n  x  k  matrix, 
and  V  is  an  k  x  n  matrix.  If  the  rank-k  update 

(2.69)  (A  +  UV*) 

of  the  matrix  A  is  invertible,  then  its  inverse  is 

(. A  +  UV*)-1  =  A-1  -  A~lU(I  +  VtA~1U)~1VtA~1. 


In  this  paper,  given  an  m  x  n  real  matrix  A,  we  will  utilize  the  Sherman- Morrison- 
Woodbury  formula  to  perform  a  specific  type  of  update  to  the  inverse  of  the  m  x  m 
matrix  AAL  In  particular,  we  wish  to  update  the  inverse  of  AAl  in  order  to  form  the 
inverse  of  the  matrix  BBt ,  where  B  is  obtained  from  A  by  deleting  its  jth  column. 
That  is, 

B  =  A  —  uv*, 

where  u  is  the  m  x  1  vector  which  is  the  jth  column  of  the  matrix  A  and  v  is  the 
n  x  1  vector  with  entries 


(2.70) 


if  i  =  j 
otherwise. 


A  simple  calculation  shows  that 

(2.71)  (A  -  uv *)  (A  -  ui^y  =  AAf  -  uvtAt  -  Aulv  +  uifvu* 

(2.72)  =  AAt  -  uu\ 


In  other  words,  if  we  form  the  matrix  B  by  deleting  the  jth  column  of  the  matrix  A, 
then  BBl  can  be  formed  from  AAf  by  a  rank-1  update.  The  Sherman-Morrison- 
Woodbury  formula  then  implies  that  the  inverse  of  BBf  can  be  computed  from 
(AA4)-1  via  a  rank-1  update. 
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3.  Numerical  apparatus 

3.1.  Nested  Legendre  discretizations  of  sequences  of  functions.  In  this  paper, 
we  are  confronted  with  sequences  (f>i, . . . ,  (pm  of  functions  defined  on  the  interval  [a,  b] 
for  which  we  wish  to  construct  a  generalized  Gaussian  quadrature  rule.  The  collection 
of  functions  (f>j  has  the  following  properties: 

•  Each  function  <pj  is  integrablc  on  [a,  6]  and  analytic  except  at  a  small  number 
of  points, 

•  The  total  number  functions  is  quite  large  (e.g.,  m  =  10,000), 

•  The  rank  of  the  set  0i, . . . ,  (pm  is  low  (e.g.,  40)  to  high  precision. 

In  order  to  construct  an  efficient  quadrature  rule,  we  will  take  advantage  of  rank 
deficiency  to  “compress”  the  sequence  <pi, . . .  ,(pm.  This  means,  more  specifically, 
that  we  first  form  a  set  Ui, . . . ,  of  orthonormal  functions  on  [a,b]  such  that  each 
pi  can  be  approximated  by  a  sum  of  the  form 

k 

(3.1)  (f>i(x)  ~  ^  a 

3  =  1 

and  each  Uj  is  defined  by  a  sum 

m 

(3.2)  Uj{x)  =  y 

i= 1 

We  then  build  a  generalized  Gaussian  quadrature  rule  for  the  functions  u\, . . .  ,Uk- 
In  the  course  of  constructing  that  quadrature,  it  will  be  necessary  to  evaluate  the 
functions  Uj  at  arbitrary  points  on  the  interval  [a,  b\.  If  the  sums  (3.2)  involve  a  large 
number  of  terms,  or  if  the  evaluation  of  the  pi  is  expensive,  then  it  is  impractical  to 
use  formula  (3.2)  to  evaluate  the  iq .  It  will  therefore  be  necessary  to  represent  the 
function  Uj  in  a  manner  which  allows  for  their  rapid  evaluation. 

An  obvious  alternative  to  evaluating  sums  of  the  form  (3.2)  directly  is  to  rep¬ 
resent  the  functions  Uj  via  polynomial  interpolation.  Let  xi,...,xn  be  a  mesh  of 
interpolation  points  on  the  interval  [a,  b]  and  suppose  that  the  Lagrange  polynomials 
interpolating  u\, ...  ,Uk  at  these  mesh  points  represent  the  functions  Uj  to  a  given 
precision.  Then,  clearly,  the  Lagrange  polynomial  interpolating  the  function 

n 

(3.3)  Vi 

3= 1 

at  xi, ...  ,xn  approximates  <pi  with  controllable  precision.  As  was  discussed  in  Sub¬ 
section  2.3,  if  Gaussian  nodes  are  used  as  interpolating  points,  then  the  resulting 
procedure  is  numerically  stable.  However,  when  the  functions  <pj  are  not  smooth, 
polynomial  interpolation  becomes  inefficient  and  can  fail  completely  for  sufficiently 
singular  functions. 

In  such  cases,  where  the  functions  to  be  interpolated  are  singular,  it  is  customary  to 
use  an  adaptive  interpolation  scheme  instead.  That  is,  the  interval  [a,  b]  is  subdivided 
into  a  collection  of  subintervals  such  that  each  function  <pj  is  accurately  interpolated 
by  a  low  order  polynomial  on  each  subinterval. 

In  this  subsection,  we  describe  a  numerical  procedure  for  the  approximation  of 
a  sequence  of  functions  via  nested  Gauss-Legendre  polynomial  interpolation.  The 
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procedure  is  introduced  in  [20],  but  we  reproduce  it  here  since  it  is  an  integral  part 
of  the  algorithm.  The  input  to  this  procedure  is  a  collection  fa, ,  <j)m  of  functions 
integrable  on  [a,  b],  a  precision  e  >  0,  and  a  reasonably  large  integer  k  which  controls 
the  number  of  interpolation  nodes  used  on  each  subinterval  (for  the  computations  in 
this  paper,  we  used  k  =  30).  The  algorithm  proceeds  in  three  stages. 

Stage  1.  In  the  first  stage,  the  following  procedure  is  used  to  discretize  each  of  the 
fa  separately.  That  is,  the  following  sequence  of  steps  is  repeated  for  i  —  1, . . . ,  m: 

1.  Construct  the  2k  Legendre  nodes  x\, ... ,  X2k  on  the  interval  [a,  b], 

2.  Let  P(x)  denote  the  Lagrange  polynomial  of  order  2k  —  1  interpolating  <f>i  at 
the  mesh  points  X\, . . . ,  x2k-  Construct  the  coefficients  ai , . . , ,  a2k  of  P(x)  in  the 
expansion 

2fc— 1 

P(X)  =  a3LAX) 
j= 0 

where  Lj{x)  is  the  jth  order  Legendre  polynomial. 

3.  If  the  inequality 

2k 

(3.4)  5]  M2<e 

j=k+l 

is  satisfied,  then  we  conclude  that  the  order  k  Legendre  expansion  for  fa  on  the 
interval  [a,  b]  is  sufficient.  Otherwise,  we  split  the  interval  [a,  b]  into  two  subintervals, 
[a,  (a  +  b)/2]  and  [(a  +  b)/ 2,  6)],  and  repeat  the  procedure  recursively  for  each  of  the 
subintervals. 

Stage  2.  Store  the  endpoints  of  each  subinterval  constructed  in  Stage  1  in  an  array. 
Sort  the  array  and  eliminate  multiple  elements.  The  resulting  array  of  points  on  the 
interval  [a,  b]  is  the  array  of  endpoints  of  the  subintervals  of  the  final  subdivision. 

Stage  3.  Construct  the  k  point  Legendre  discretization  on  each  the  subintervals 
obtained  in  Stage  2  for  each  of  the  functions  fa(x). 

Remark  3.1.  The  scheme  of  this  subsection  is  a  reasonably  reliable  mechanism  for 
the  discretization  of  sets  of  functions  with  singularities.  The  stopping  condition  (3.4) 
is  analogous  to  that  usually  used  to  terminate  an  adaptive  quadrature  procedure.  Just 
as  any  such  quadrature  procedure  can  fail  for  carefully  designed  counterexamples,  so 
too  can  the  procedure  of  this  section  fail  under  certain  circumstances.  The  problem, 
however,  is  not  encountered  in  the  examples  of  this  paper  and  whenever  the  authors 
have  encountered  it  in  practice,  it  has  been  easy  to  rectify. 

3.2.  Compression  of  a  finite  sequence  of  functions.  This  subsection  describes  a 
numerical  procedure  for  compressing  a  finite  sequence  of  functions  by  approximating 
either  the  singular  value  decomposition  or  a  QR  decomposition  of  the  sequence. 

We  begin  with  two  definitions,  the  second  of  which  is  adapted  from  [20]  and  [3] .  In 
what  follows,  PC([a,b\ )  refers  to  the  vector  space  of  piecewise  continuous  functions 
on  the  interval  [a,  b]. 

DEFINITION  3.1.  A  k-point  linear  interpolation  scheme  on  the  interval  [a,  b]  is 
a  collection  of  linearly  independent  functions  fa,  ...fa  in  PC([a,b ])  and  a  set  of 
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distinct  nodes  X\,...,Xk  in  [a,b]  together  with  a  linear  mapping  T  :  PC([a,b ])  — > 
span{4> i,  •  •  • ,  4>k}  such  that 

(3.5)  Tf(Xj)  =  f{xj ) 

for  all  j  =  1, . . . ,  k. 

We  call  the  x±, ...  ,Xk  interpolation  nodes,  the  mapping  T  the  interpolation  map¬ 
ping,  and  the  <pj  interpolating  functions.  We  also  say  that  the  coefficients  aq , ,atk 
ofTf(x)  with  respect  to  the  basis  {4>i, . . .  ,4>k}  are  the  interpolation  coefficients  for 
the  function  f . 

Associated  with  every  fc-point  linear  interpolation  scheme  is  an  invertible  linear 
transformation  — >  Mfc  which  takes  the  values  f(xi),...,f(xk)  of  a  function  / 
at  the  interpolation  nodes  to  the  k  interpolation  coefficients  for  the  function.  The 
image  Tf  of  /  under  the  interpolation  mapping  is,  of  course,  determined  by  these 
interpolation  coefficients.  We  will  often  refer  to  Tf  as  the  function  defined  by  the 
interpolation  scheme  and  the  values  f(x i), . . . ,  f(x]f). 

Remark  3.2.  It  is  generally  possible  to  extend  the  definition  of  interpolation  scheme 
to  encompass  interpolating  functions  f>  not  in  PC([a,b}),  as  well  as  interpolation 
mappings  T  defined  on  larger  spaces.  We  must  keep  in  mind,  however,  that  equation 

(3.5)  requires  that  both  T f{x)  and  f(x)  be  defined  pointwise. 

DEFINITION  3.2.  We  say  that  the  the  combination  of  a  k-point  linear  inter¬ 
polation  scheme  on  [a,  b]  with  nodes  y\, . . . ,  e  [a,  b]  and  interpolating  functions 
4> i, . . .  ,f>k,  and  a  k-point  quadrature  rule  with  nodes  aq, . . .  ,Xk  G  [a,  b]  and  weights 
Wi, ...  ,Wk  preserves  inner  products  if 

•  The  nodes  x\, ...  ,Xk  of  the  quadrature  rule  coincide  with  the  nodes  y\, ...  ,yk 
of  the  interpolation  scheme,  and 

•  The  quadrature  rule  is  exact  for  all  products  4>i(x)4>j{x )  of  the  interpolating 
functions. 

Remark  3.3.  The  second  condition  in  definition  3.2  together  with  the  assumption  of 
the  linearity  of  the  interpolation  scheme  ensures  that  for  any  two  piecewise  continuous 
functions  f(x)  and  g(x)  the  quadrature  rule  is  exact  for  the  integral 

(3.6)  /  T f(x)Tg(x)dx. 

J  a 

The  following  are  examples  of  quadrature  and  interpolation  schemes  which  preserve 
inner  products: 

Example  3.1.  The  combination  of  a  classical  Gaussian  quadrature  and  Lagrange 
interpolation  at  the  same  Gaussian  nodes  preserves  inner  products,  since  polynomials 
interpolation  on  n  nodes  produces  an  interpolating  polynomial  of  order  n  —  1  and 
the  product  of  any  two  such  polynomials  is  exactly  integrated  by  an  n  point  Gaussian 
quadrature. 

Example  3.2.  Nested  Gaussian  interpolation  —  like  that  of  the  preceding  section  - 
coupled  with  corresponding  nested  Gaussian  quadrature  formulas  also  preserve  inner 
products,  since  on  each  subinterval  Gaussian  interpolation  is  coupled  with  a  Gaussian 
quadrature  formula. 

Note  that  the  interpolating  functions  in  this  example  are  not  continuous,  but  rather 
piecewise  continuous. 


18 


The  following  theorem,  which  is  a  reformulation  of  Theorem  4.5  in  [20],  will  be 
used  in  approximating  the  singular  value  decomposition  (or  a  QR  decomposition)  of 
a  sequence  of  functions. 

THEOREM  3.1.  Suppose  that  the  combination  of  a  quadrature  rule  with  nodes 
Xi, . . .  ,xn  and  weights  uq, . . . ,  wn,  and  a  linear  interpolation  scheme  with  interpo¬ 
lating  functions  ifi, ...  ,ifn  and  interpolation  mapping  T  preserves  inner  products  on 
[a,  b\.  Suppose  further  that  4>i, . . . ,  0m  is  a  collection  of  piecewise  continuous  func¬ 
tions  on  [a,  b\,  U  —  ( Uij )  is  an  nx  k  matrix  with  orthonormal  columns,  mid  R  =  (rtj) 
is  a  k  x  m  matrix  such  that 

k 

(3.7)  (j)j{Xi)y/wl  =  y,  UuXlj 

1=1 


for  all  j  =  1, . . .  ,m  and  i  —  1, ...  ,n.  If  the  functions  Uj(x )  are  defined  for  all 
j  =  1 , ...  ,k  via  the  interpolation  scheme  and  the  values 


(3.8) 


Uj(Xi) 


then: 

1.  The  functions  Uj(x )  are  orthonormal;  that  is, 


Ui(x)uj(x)dx  =  Sij 

for  all  i,j  =  1, . . . ,  k. 

2.  The  sequence  of  functions  <f>\, ... ,  <f>m  defined  by  the  formula 

k 

<j>j{x)  =  ui(x)rij 

i=  1 

is  identical  to  the  sequence  of  functions  produced  by  sampling  the  functions 
<f>i, ... ,  (fim  at  the  points  {aq}  and  then  interpolating  with  the  interpolation  scheme 
on  [a,  b] . 

The  algorithm  described  below  uses  as  input  a  sequence  of  functions  <t>u  . .  -  ,  <pm  '■ 
[a,  b]  — >  M  and  a  quadrature  and  interpolation  scheme  which  preserves  inner  products. 
The  nodes  and  weights  of  the  quadrature  will  be  denoted  by  x\r.  . .  ,  xn  and  w\, ... ,  wn, 
respectively. 

The  output  of  the  algorithm  depends  on  whether  the  SVD  or  a  QR  decomposi¬ 
tion  is  used  for  compression.  In  either  case,  the  output  is  a  sequence  u\, ...  ,Uk  of 
orthonormal  vectors  and  a  sequence  a i, . . .  ,<7*.  of  positive  real  numbers.  When  the 
SVD  is  used,  the  orthonormal  vectors  u±, . . . ,  Uf.  approximate  the  singular  vectors  for 
the  sequence  4>i,  ■  ■  • ,  (ftm  and  the  ai, ...  ,ak  approximate  the  corresponding  singular 
values.  In  the  case  of  the  QR  decomposition,  the  orthonormal  vectors  U\, ...  ,uk  ap¬ 
proximate  a  collection  of  QR  vectors  for  4>i,  •  •  • ,  4>m  and  the  oq  are  the  corresponding 
normalizing  factors. 

Description  of  the  algorithm. 

1.  Construct  the  n  x  m  matrix  A  with  entries 

■A-ij  0j  (■-*'/. )  \J  a')  • 
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2.  Compute  either  the  SVD  of  the  matrix  A  or  a  QR  decomposition  for  the  matrix 
A.  In  the  first  case,  produce  the  factorization 


A  =  UEV*, 

where  U  =  (■ Uij )  is  an  n  X  k  matrix  with  orthonormal  columns,  V  =  (vtj)  is  a  n  mx  k 
matrix  with  orthonormal  columns,  and  S  is  a  diagonal  matrix  whose  jth  diagonal 
entry  is  ar  In  the  second  case,  produce  the  factorization 

ATI  =  UR, 


where  U  =  (ul3)  is  an  n  x  k  matrix  with  orthonormal  columns  and  R  is  an  k  x  m 
trapezoidal  matrix  with  diagonal  entries  <J\, . . . ,ak. 


3.  Construct  the  n  x  k  values  Uj(xf)  defined  by  the  formula 


(3.9) 


u. 


Uj{Xi)  = 


V 


4.  For  any  desired  point  x  e  [o,  6],  evaluate  the  functions  rq  :  [ a,b ]  — >  R  using  the 
interpolation  scheme  on  [a,  b]. 

Remark  3.4.  Theorem  3.1  ensures  that  the  accuracy  of  the  approximation  produced 
by  the  algorithm  of  this  section  depends  primarily  on  the  accuracy  of  the  underlying 
interpolation  scheme. 


3.3.  Construction  of  Chebyshev  quadratures.  In  this  subsection,  we  describe  a 
numerical  algorithm  for  the  construction  of  a  Chebyshev  quadrature  for  a  finite  se¬ 
quence  of  functions.  Given  a  pre-existing  n-point  quadrature  formula  x\,  x-i, . . . ,  xn ,  w\,  w\, . . . ,  wn, 
where  n  >  k,  which  exactly  integrates  the  U\, . . . ,  Uk ,  it  is  straightforward  to  construct 
a  Chebyshev  quadrature  for  u±, . . .  ,Uk .  Because  the  pre-existing  quadrature  rule  ex¬ 
actly  integrates  the  input  functions  U\, . . .  ,Uk,  the  matrix  equation 


/  Ui(x\)  Ui(x2)  •••  Ui(xn)  \ 

f  w  1  \ 

/  ri  \ 

U2(x  i)  u2{x2)  •••  u2{xn) 

w2 

— 

r2 

\  uk(x i)  uk(x2 )  •  •  •  uk(xn )  / 

\  Wn  J 

\  rk  ) 

where  rl  is  dehned  by 


for  all  i  =  1, . . . ,  k,  is  satisfied, 
such  that 


fb 

rt  =  /  Uj{x)dx 
J  a 

By  Lemma  2.5,  there  exist  i i, 


in  and  w\, . .  .vdk 


(3.11) 


/  ui(z*i) 

u2(xi2 

life 

Ul  (xin)  \ 
MXin ) 

Uk{xin)  / 


li\ 


wk 

0 


/  n  \ 

r2 

\  n  ) 


\  0  ) 
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and 

k  n 

(3.12)  ^  \wj\  <  k  \wj\. 

1=1  1=1 

In  other  words,  there  is  by  Lemma  2.5  a  /c-point  quadrature  formula  xtl , ... ,  xik,  Wi, . . . 
for  the  k  functions  ui,...,uk.  Moreover,  the  inequality  (3.12)  ensures  that  the  k- 
point  quadrature  formula  is  numerically  stable  assuming  the  weights  of  the  original 
quadrature  are  reasonably  small. 

The  following  algorithm  is  a  computational  procedure  for  constructing  such  a  quad¬ 
rature  via  the  modified  Gram-Schmidt  algorithm  with  double  orthogonalization,  or 
one  of  its  variants.  This  procedure  does  not  realize  the  theoretical  bound  (2.5)  guar¬ 
anteed  by  Lemma  2.5.  In  practice,  however,  it  still  leads  to  the  construction  of  stable 
Chebyshev  quadrature  formulas  for  the  input  functions;  see  Subsection  2.8  for  a  more 
discussion  of  the  numerical  stability  of  the  modified  Gram-Schmidt  algorithm  and  its 
variants. 

The  algorithm  uses  as  input  a  sequence  Ui,...,uk  of  functions  and  a  pre-existing 
quadrature  x±, ... ,  xn,  w i, ... ,  wn,  with  n  >  k,  which  exactly  integrates  the  functions 
ui,  ■  ■  ■ ,  Uk,  and  produces  as  output  a  /e-point  quadrature  formula  consisting  of  k  nodes 
£i, . . . ,  Xk  G  {x\, . . . ,  xn}  and  k  weights  w\, . . . ,  Wk- 

Description  of  the  algorithm. 

1.  For  the  vector  r  G  whose  ith  entry  is  the  sum 

n 

r*  =  ^Uiix^Wj, 

3= i 

which  is  the  value  of  the  integral 

Ui{x)dx. 

2.  Form  the  k  x  n  matrix  B  with  entries 

Bij  =  Ui(Xj)y/w]. 

3.  Lise  the  pivoted  Gram-Schmidt  algorithm  with  double  reorthogonalization  to  select 
a  set  Btl , ,  Bik  of  spanning  columns  for  the  matrix  B  and  form  the  factorization 

£)• 

where  Rn  is  a  k  x  k  upper  triangular  matrix,  Q  is  a  k  x  k  orthogonal  matrix, 

(  Bi i  BV2  ■  ■  ■  Bik  )  =  QRu, 

and  || i?22 1| 2  is  small. 

4.  Lise  back  substitution  to  construct  a  solution  z  G  to  the  k  x  k  system  of 
equations 

Rnz  =  Q*r, 

which  is,  of  course,  a  least  squares  solution  of  the  system 

(  Bh  Bh  ■  ■  ■  Bik  )  z  =  r. 


,wk 
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5.  Form  the  new  h-point  quadrature  £\, . . . ,  Xk ,  hq, , . . ,  hq  by  letting 

x,-  =  x*.  and  Wj  =  Zj^/uJ] 

for  all  j  —  1, . . . ,  k. 

Remark  3.5.  The  columns  of  the  matrix  B  are  scaled  by  the  square  roots  of  the 
quadrature  weights  so  that  the  l 2  norms  of  columns  of  B  computed  in  Step  3  by 
the  Gram-Schmidt  orthonormalization  procedure  are  proportional  to  the  quadrature 
weight  for  the  corresponding  column. 

4.  Numerical  algorithm 

This  section  describes  a  numerical  algorithm  for  the  construction  of  the  nodes 
and  weights  of  an  approximate  quadrature  rule  for  a  sequence  of  functions.  The 
algorithm’s  input  is  a  sequence  of  functions 

(4.1)  0i,  ■  •  -Am  :  [a,  b]  -*■  M 

and  two  accuracies,  €diSC  and  equad ■  The  first,  e<jjSC,  controls  the  accuracy  of  the  scheme 
used  to  discretize  and  compress  the  input  functions  0i, . . . ,  0m,  and  the  second,  equad, 
is  the  desired  accuracy  for  the  quadrature  rule.  The  algorithm’s  output  is  a  quadra¬ 
ture  rule  consisting  of  a  set  of  nodes  x\ , ,xi  and  a  set  of  weights  w\ , ...  ,wi  such 
that 

/b  l 

4>i{x)dx  ~  y^^jjx^Wj 

3= 1 

for  all  i  =  l, ...  pm.  The  integer  l  depends  on  the  numerical  rank  of  the  input  set 
0i,  •  •  • ,  0m  and  both  edisc  and  equad. 

The  algorithm  proceeds  in  three  stages.  In  the  first  stage,  the  numerical  techniques 
of  the  preceding  section  are  used  to  discretize  and  compress  the  functions  0i, . . . ,  0m. 
In  the  second,  a  suboptimal  quadrature  rule  is  obtained  for  the  compressed  collection 
of  functions.  Finally,  in  the  third  phase,  an  optimization  procedure  is  used  to  reduce 
the  number  of  points  needed  by  the  quadrature. 

Stage  1:  Discretization  and  compression. 

In  this  stage  the  following  sequence  of  steps  is  performed  to  discretize  and  compress 
the  input  functions. 

1.  Use  the  technique  of  Section  3.1  to  discretize  the  functions  0!, . . . ,  0m,  so  that  they 
are  all  represented  to  the  precision  e^sc.  Let  0i(x), . . . ,  0m(x)  denote  the  resulting 
discretizations.  Also,  let  x\, . . . ,  xn  denote  the  discretization  nodes  and  w\ , ,wn 
the  weights  of  the  associated  quadrature. 

2.  Apply  the  procedure  of  Subsection  3.2  to  compress  the  sequence  0i,...,0n  via 
either  a  QR  decomposition  or  the  singular  value  decomposition.  Denote  the  re¬ 
sulting  orthonormal  functions  by  Ui, ...  ,up  and  the  associated  singular  values  or 
normalization  factors  by  Ai  >  A2  >  . . .  >  Xp  >  0. 

3.  Discard  any  of  the  functions  u±, ...  ,up  corresponding  to  a  singular  value  (or  nor¬ 
malization  factor)  Xj  <  equad- 

4.  Denote  the  functions  obtained  in  this  manner  by  uii...Ju/e  and  the  associated 
singular  values  (normalization  factors)  by  Ai  >  A2  >  . . .  >  A*,  >  0. 
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Remark  4.1.  It  is  important  that  equad  be  somewhat  larger  than  e^sc  (for  the  exam¬ 
ples  of  this  paper,  equad/£disc  >  100,).  Because  the  functions  (p i,  •  •  • ,  <pm  are  discretized 
to  precision  ediSC,  singular  vectors  Uj  corresponding  to  singular  values  Xj  comparable 
to  edisc  contain  little  information  about  the  functions  (p\, . . .  ,(pm-  This  means  that 
constructing  a  quadrature  formula  which  faithfully  integrates  such  singular  functions 
is  pointless,  and  moreover,  since  the  Uj  corresponding  to  singular  values  A,  compa¬ 
rable  to  esvd  are  not  necessarily  smooth,  even  when  the  the  functions  (pi, ... ,  <pm  are 
very  smooth,  it  can  cause  difficulties  with  the  the  algorithm  described  in  Step  3  below. 

Stage  2:  Construction  of  a  /c-point  quadrature  rule. 

We  now  apply  the  procedure  of  Subsection  3.3  to  construct  a  /c-point  quadrature 
formula  for  (pi, ... ,  (f)m.  We  use  as  inputs  to  that  procedure  the  functions  U\, ...  ,Uk 
and  the  nested  Gaussian  quadrature  aq, . . .  ,xn,Wi , . . .  ,wn  constructed  in  Stage  1. 
Note  that  since  the  ui,...,Uk  are  defined  via  the  nested  Gaussian  interpolation 
scheme,  the  quadrature  aq, . . . ,  xn ,  uq, . . . ,  wn  is  exact  for  the  functions  Ui, . . . ,  Uk . 

Denote  the  resulting  k  quadrature  nodes  by  aq, . . . ,  £k  and  the  k  quadrature  weights 
by  wi, . . . ,  Wk-  The  weights  of  the  nested  Gaussian  quadrature  scheme  satisfy  the 
bound 

n 

(4.3)  <  (6-  a), 

3= 1 

since  w3  >  0  for  all  j  =  1 , ...  ,n  and  the  quadrature  rule  is  exact  for  the  function 
f(x)  =  1.  It  follows  from  the  discussion  in  Subsection  3.3  that  the  new  /c-point 
quadrature  formula  is  numerically  stable 

Because  the  discretizations  (p  i, . . .  ,<f>m  are  well  approximated  by  functions  in  the 
span  of  the  u±, . . .  ,Uk,  a  quadrature  rule  exact  for  the  functions  Ui, . . .  ,Uk  will  ap¬ 
proximately  integrate  the  functions  (pi, ..  . . ,  (j)m  (and,  of  course,  it  follows  that  such 
a  quadrature  rule  will  also  serve  for  the  functions  <pj{x ),  assuming  the  discretiza¬ 
tions  constructed  in  Step  1  are  sufficiently  accurate).  So  the  new  /c-point  quadra¬ 
ture  rule  xi, . . .  ,Xk,wi, . . .  ,Wk  is  an  approximate  quadrature  for  the  input  functions 
01,  .  .  .  ,  (pm  ■ 

Stage  3:  Point-by-point  reduction  of  the  quadrature  rule. 

In  this  stage,  beginning  with  an  n-point  quadrature  rule  aq, . . . ,  xn,Wi, . . . ,  wn,  we 
repeatedly  apply  the  following  sequence  of  steps,  which  attempt  to  reduce  an  n-point 
quadrature  rule  for  u\, . . .  ,Uk  to  an  (n  —  l)-point  quadrature  rule  for  u\, . . . ,  Uk ■ 

To  perform  these  calculations,  we  use  the  observation  of  Subsection  2.2;  namely, 
that  the  nodes  and  weights  of  a  quadrature  rule  exact  for  a  sequence  of  functions 
satisfy  a  system  of  nonlinear  equations.  Given  one  of  the  n  quadrature  nodes  x3,  we 
can  use  this  fact  to  compute  an  (n  —  1)  point  quadrature  for  the  functions  Ui, ...  ,Uk 
via  the  damped  Gauss-Newton  algorithm  described  in  Subsection  2.4:  we  simply  form 
the  system  of  k  nonlinear  equations  in  n  —  1  unknowns  described  in  Subsection  2.2 
and  use  as  an  initial  guess  for  the  Gauss-Newton  iterations  the  current  quadrature 
nodes  and  weights,  excluding  the  chosen  point  Xj  and  its  corresponding  weight  Wj. 

The  only  problem  with  this  approach  is  the  selection  of  the  quadrature  node  Xj 
to  eliminate.  There  are  n  possible  nonlinear  systems,  each  corresponding  to  one  of 
the  n  points  which  can  be  omitted.  It  is  computationally  expensive  to  search  for  a 
sufficiently  accurate  (n  —  l)-point  quadrature  rule  by  solving  each  of  the  n  resulting 
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nonlinear  systems  of  equations  via  the  damped  Gauss-Newton  method.  That  dif¬ 
ficulty  is  surmounted  via  the  procedure  described  in  Step  1,  wherein  the  direction 
of  the  step  for  the  first  iteration  of  the  damped  Gauss-Newton  method  is  computed 
for  each  of  the  n  possible  nonlinear  systems.  The  procedure  is  expedited  by  using 
the  Sherman-Morrison- Woodbury  update  formula,  which  leads  to  a  computationally 
feasible  scheme,  and  the  results  are  used  to  determine  the  order  in  which  to  remove 
quadrature  nodes. 

Step  1:  Rank  the  remaining  nodes. 

1.  For  each  i  =  1, . . . ,  k,  compute  the  integrals 


u.i(x)dx 

using  the  original  nested  Legendre  quadrature  rule  formed  in  Step 
vector 


/  n  \ 

G 

\  rk  ) 


1. 


Form  the 


2.  Form  the  Jacobian  matrix 

/  ui(xi)wi  u'1(x2)w2 


J  = 


u'2(xi)wi  u'2{x  2)w2 


u[(xn)wn  Ui{x\)  Ui(x2) 
u2(xn)wn  u2(x  i)  u2{x2) 


Ul(xn)  \ 
u2{xn) 


\  u'k(x i)wi  uk(x2)w2  . . .  u'k{xn)wn  uk{x i)  uk{x 2)  . . .  un{xn)  ) 
for  the  nonlinear  system  (2.4)  in  Subsection  (2.2),  and  compute  the  inverse 

A  =  (jfy1. 


of  the  product  JJ*. 

3.  For  each  node  xk,  hrst  use  the  Sherman- Morrison- Woodbury  formula  (see  Subsec¬ 
tion  2.9)  to  form  the  matrix 

Ak  =  (JkJkT1 


via  two  rank-1  updates  to  A,  where  the  matrix  Jk  is  obtained  from  J  by  deleting 
its  kth  and  (k  +  n)th  columns.  That  is,  Jk  is  obtained  from  J  by  deleting  the 
contributions  from  the  node  xk  and  its  corresponding  weight  wk. 

Next,  compute  the  damped  Gauss-Newton  step  direction  Axk  for  the  nonlinear 
system  obtained  by  omitting  the  point  xk,  i.e.,  find  a  solution  Axk  to  the  least 
squares  problem 

argmin  ||  Jkx  —  r\\2, 

X 

where  r  is  the  vector  of  definite  integrals  computed  above.  The  solution  to  the 
minimization  problem  is  found  by  solving  the  normal  equations 

Axk  =  {JkJlY1  J{r  =  AkJlkr. 
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4.  For  each  k  =  1, ...  ,  n,  let  rjp.  denote  the  l2  norm  of  the  solution  vector  Aaq ■  We 
will  refer  to  the  value  rjk  as  the  significance  of  the  node  aq.. 

5.  Renumber  the  nodes  and  weights  of  the  quadrature  so  that  {aq, . . . ,  xn }  are  ar¬ 
ranged  in  order  of  increasing  rjj. 

Step  2:  First  pass  through  the  nodes. 

For  each  j  =  1, ...  ,n,  perform  the  following  sequence  of  steps: 

1.  Form  an  initial  guess  aq, . ... . ,  Xj, . . . ,  xn  and  uq , . . .  ,wj, . . .  ,wn  for  the  damped 
Gauss-Newton  method,  which  excludes  the  node  x3  and  its  corresponding  weight 

Wj. 

2.  Perform  a  small  number  (e.g.,  4)  of  damped  Gauss-Newton  iterations  to  form  a 
(n— l)-point  quadrature  rule  aq, . . . , Xn-i,  uq, . . . ,  Wn- i  for  the  functions  u±, . . . ,  Uk- 

3.  Measure  the  approximation  error 

k  n—  1 

g  =  5Z  - ri 

i=l  j= 1 

for  this  new  quadrature  rule. 

4.  If  the  error  e3  for  the  quadrature  rule  is  sufficiently  small  (i.e.,  e3  <  equad),  then  we 
accept  this  (n  —  l)-point  quadrature  rule  and  go  to  Step  4. 

If  this  sequence  of  steps  completes  without  finding  a  quadrature  rule  with  acceptable 
accuracy,  then  we  renumber  the  points  x\, ...  ,xn  and  weights  uq , . . . ,  wn  so  that 

G  <  ^2  <  •  •  •  <  €n, 

and  move  onto  Step  3. 

Step  3:  Second  pass  through  the  nodes. 

We  arrive  at  this  stage  only  if  we  were  unable  to  find  a  satisfactory  (n  —  l)-point 
quadrature  rule  by  taking  a  small  number  of  damped  Gauss-Newton  steps.  For  each 
j  =  1 , . . . ,  n  we  perform  the  following  sequence  of  steps: 

1.  Form  an  initial  guess  aq, . . . ,  Xj, _ ,  xn  and  «q, . . . ,  Wj, . . . ,  wn  for  the  clamped 

Gauss-Newton  method,  which  excludes  the  node  Xj  and  its  corresponding  weight 
w3. 

2.  Use  the  damped  Gauss-Newton  algorithm  to  form  an  (n  —  l)-point  quadrature 
rule  xi, . . . ,  x^- 1,  Wi, . . . ,  w^~ i  for  the  functions  Ui, . . . ,  Uk-  In  this  step,  the 
limit  m  on  the  number  of  iterations  should  be  large  (for  the  examples  of  this 
paper,  m  =  30). 

3.  Measure  the  approximation  error 

k 


2=1 

for  this  new  quadrature  rule. 

4.  If  the  error  e3  for  the  quadrature  rule  is  sufficiently  small  (i.e.,  e3  <  equad),  then 
we  accept  this  (n  —  l)-point  quadrature  rule  and  go  to  Step  4. 

Step  4:  Form  the  (n  —  l)-point  quadrature. 

If  an  (n  —  1)  point  quadrature  rule  with  sufficient  precision  has  been  found,  then 
we  accept  this  rule  and  repeat  the  procedure  of  this  stage  for  the  newly  formed 
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Xi 

Wi 

Xi 

Wi 

0.2768757118897219E-21 

0.6728999118496260E-16 

0.2913313729367070E-12 

0.1363525818876240E-09 

0. 1466530079834641E-07 
0.5526176228089556E-06 
0.9556438762298892E-05 
0.9137352164364451E-04 
0.5528197186769772E-03 
0.2340613778507800E-02 
0.7487402952102819E-02 
0.1919273618673196E-01 
0.4125089438263980E-01 
0.7699431299223189E-01 

0. 1282821075952860E+00 

0. 1950313100451387E+00 
0.2754007895160087E+00 
0.3663421175354336E+00 
0.4641496390531796E+00 
0.5648167175087882E+00 
0.6641918632800308E+00 
0.7580165934350541E+00 
0.8419459185715879E+00 
0.9116491627096593E+00 
0.963071 1719694888E+00 
0.9928824356753315E+00 

0.4385873207003101E- 
0.6670167705957102E- 
0.2065563299999408E- 
0.7260286603891295E- 
0.6011573165242795E- 
0.1770124736858168E- 
0.2415867272187855E- 
0. 1836327943155820E- 
0.8884287232428108E- 
0.3023731195124369E- 
0.7813101428951048E- 
0.1624899441619337E- 
0.2844094097856474E- 
0.4336382757029246E- 
0.5919402249047975E- 
0.7398735448703726E- 
0.8623391772361065E- 
0.9502522703162401E- 
0.9991966301938738E- 
0.1007221934294668E- 
0.9731824512000341E- 
0.8960557633305366E- 
0.7752702387328704E- 
0.6119598511874399E- 
0.4108284432044861E- 

0. 1819930030488772E- 

200.: 

15)4 

in.: 

osd.: 

on.: 

05). ‘ 

04).; 

03).; 

03).; 

02)4 

02)4 

on.: 

on.^ 

on.) 

Offl.l 

00.1 

00.2 

00.3 

00.4 

005 

00.6 

00.7 

00.8 

00.9 

00.9 

00.9 

056470082726196E-21 
448276636263963E-16 
844497556672273E-12 
004613396825143E-09 
199010226476309E-07 
834665505524823E-06 
721446282278398E-05 
558255329196107E-04 
262332881 183370E-03 
251655341009627E-02 
255657724649352E-02 
869906517557028E-01 
035483211529235E-01 
555760069233022E-01 
261876307279175E+00 
)22011617225488E+00 
718197580245527E+00 
320819036055107E+00 
593842035808489E+00 
598231960985849E+00 
593340785103139E+00 
537053806233217E+00 
585666970462145E+00 
)94519100536209E+00 
320432996182802E+00 
)26707956133087E+00 

0. 1711096457804046E-20 
0.3503308273337532E-15 
0.1340813788691161E-11 
0.5475770404764035E-09 
0.5015698242641016E-07 

0. 1574374049854322E-05 
0.2233377157238088E-04 

0. 1737239626794526E-03 
0.8524600446470783E-03 
0.2928241357781042E-02 

0.7616225404410970E-02 

0. 1591836338093090E-01 
0.2796728481689857E-01 
0.4275972822079389E-01 
0.5848899820837419E-01 
0.7323146146527266E-01 
0.8550269331832136E-01 
0.9441537793180152E-01 
0.9953654606755140E-01 

0. 1006639600889222E+00 
0.9765633871558868E-01 
0.9035707472894443E-01 
0.7861734429197797E-01 

0. 624251 1779052747E-01 
0.4213323000112514E-01 

0.1873187135192821E-01 

TABLE  1.  Quadrature  formulas  for  the  functions  of  the  form  (5.2)  with 
a  G  [—.6, 1.0]  and  b  =  20.  The  26-point  rule  on  the  left  was  generated 
with  the  QR  variant  of  the  algorithm,  and  the  26-point  rule  on  the 
right  with  the  SVD  variant.  Both  achieve  full  double  precision  accu¬ 
racy. 


(n  —  1)  point  quadrature,  beginning  with  Step  1.  Otherwise,  we  accept  the  n-point 
quadrature  rule  which  was  the  input  to  this  stage  and  the  algorithm  terminates. 

Remark  4.2.  The  process  terminates  when  an  n-point  quadrature  rule  cannot  be 
reduced  to  an  (n  —  1  )-point  quadrature  rule  without  an  unacceptable  loss  of  precision. 

Remark  4.3.  We  would  like  to  reiterate  that  the  quadrature  rules  obtained  by  this  al¬ 
gorithm  are  not  exact  for  the  functions  4>i,  •  •  • ,  4>m,  but  rather  depend  on  the  precision 
of  the  approximation  used  and  the  pointwise  properties  of  the  (fj. 

5.  Numerical  Examples 

We  have  implemented  the  algorithm  of  this  paper  for  the  computation  of  efficient 
quadrature  rules  and  tested  it  on  a  number  of  examples.  Below  we  present  sev¬ 
eral  of  these  examples  in  order  to  demonstrate  the  performance  of  the  algorithm. 
It  was  implemented  in  Fortran  77,  compiled  with  the  Lahey-Fujitsu  FORTRAN  95, 
and  all  timings  were  measured  on  a  2.0  GHz  Intel  Core  2  Duo  processor  with  2GB 
of  RAM  (no  parallelization  was  utilized).  Where  possible,  computations  were  per¬ 
formed  in  double  precision  (FORTRAN  REAL*8)  arithmetic;  however,  in  order  to 
achieve  double  precision  accuracy  for  quadrature  rules,  it  is  necessary  to  perform  the 
computations  in  extended  precision  (FORTRAN  REAL*  16) .  Because  the  Intel  Core 
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Xi 

Wi 

Xi 

Wi 

0.7142868061585990E-22 

0.2456099244659448E-16 

0.1388771056068208E-12 

0.7923942786393112E-10 

0.9728977588843004E-08 

0.3945858602101462E-06 

0.7004999008550828E-05 

0.6647847640796673E-04 

0.3909094695073146E-03 

0.1591555788166976E-02 

0.4880565012011741E-02 

0.1201508322521625E-01 

0.2491818304126502E-01 

0.4515299730588538E-01 

0.7352628756922541E-01 

0. 1099997234334442E+00 

0. 1538773279072512E+00 
0.2040969860767333E+00 
0.2594754054419947E+00 
0.3188532789291890E+00 
0.3811550534088713E+00 
0.4453957796025275E+00 
0.5106608382341754E+00 
0.5760733409165308E+00 
0.6407559864297044E+00 
0.7037897986529295E+00 
0.7641706741669453E+00 
0.8207656867459757E+00 
0.8722752374352830E+00 
0.9172156595977674E+00 
0.9539498396044303E+00 
0.9808044766847136E+00 
0.9963048714773390E+00 

0.1160152107855531E- 
0.2509065273801540E- 
0.1016875670063657E- 
0.4349600879799666E- 
0.4086717652474466E- 
0.12827398311U961E- 
0.1775840695185165E- 
0.1322763958769828E- 
0.6144256207665815E- 
0. 1989634887947442E- 
0.4888053092794360E- 
0.9711002899742405E- 
0. 1636280699646038E- 
0.2424264956181601E- 

0.3249392541983736E- 

0.4033170405949942E- 

0.4724229693878507E- 

0.5299744130306683E- 

0.5756526955130879E- 

0.6101U6645139500E- 

0.6342859582873963E- 

0.6490097796U0366E- 

0.6548363166007094E- 

0.6519584582093280E- 

0.6401701018681496E- 

0.6188400467055524E- 

0.5868958877875049E- 

0.5428419177206927E- 

0.4848722081318521E- 

0.4111890543212653E- 

0.3206665104837302E- 

0.2139086878064731E- 

0.9450771617397615E- 

200 

L50 

L10 

)90 

)70 

)50 

)  10 

)30 

)30 

)20 

)20 

)20 

)10 

)10 

)10 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

m. 

)D. 

)D. 

)D. 

)D. 

)D. 

)29. 

9545672035004698E-22 
2782141914841101E-16 
1433442445848540E-12 
77935721 17442273E-10 
9385234534215939E-08 
379781 1590248750E-06 
6783814687534227E-05 
6495463184473910E-04 
3852594648391139E-03 
1580066444838872E-02 
4874206864113175E-02 
1205644750230908E-01 
2509445953623319E-01 
4558557083460269E-01 
7433570351228120E-01 
L112672046075921E+00 
L556251204014568E+00 
2062913928215772E+00 
2620423475006598E+00 
S216956048993508E+00 
J841661902355933E+00 
1484687328581966E+00 
U36935077033405E+00 
5789715342131883E+00 
i434353506653515E+00 
r061776301391591E+00 
i7662083501992056E+00 
^224122548819472E+00 
^735123565945287E+00 
U80534258335345E+00 
)544319708711901E+00 
)810102430752318E+00 
)963450010U5773E+00 

0.1529382502435904E-20 

0.2808039205209724E-15 

0. 103991 1929764398E-1 1 
0.4253421034031052E-09 
0.3933760334981154E-07 

0. 1235803705030886E-05 

0. 1725353179365193E-04 

0. 1298334366282352E-03 
0.6086396169406588E-03 

0. 1985699293680274E-02 
0.4907602643213011E-02 
0.9794765851446592E-02 
0.1655504783796611E-01 
0.2456292098104041E-01 

0.3292041369959096E-01 

0.4081110237571449E-01 

0.4771399178998336E-01 

0.5341223120206637E-01 

0.5789177660270919E-01 

0.6123393859233695E-01 

0.6354343717872484E-01 

0.6491062297564529E-01 

0.6539475441359519E-01 

0.6501741217619580E-01 

0.6375971643567762E-01 

0.6156055973389270E-01 

0.5831564380342947E-01 

0.5387967213251687E-01 

0.4807758967228032E-01 

0.4073548782515963E-01 

0.3174482830959985E-01 

0.2116533404178935E-01 

0.9348448278692109E-02 

TABLE  2.  Quadrature  formulas  for  the  functions  of  the  form  (5.2)  with 
a  G  [—.6, 1.0]  and  b  =  50.  The  33-point  rule  on  the  left  was  generated 
with  the  QR  variant  of  the  algorithm,  and  the  33-point  rule  on  the 
right  with  the  SVD  variant.  Both  achieve  full  double  precision  accu¬ 
racy. 


2  Dno  processor  does  not  support  extended  precision  arithmetic  in  hardware,  the 
running  times  for  these  computations  are  quite  long. 

Throughout  this  section,  we  will  denote  by  Jn(z)  the  nth  order  Bessel  function  of 
the  first  kind,  by  Yn(z)  the  nth  order  Bessel  function  of  the  second  kind,  and  by  Hn(z) 
the  nth  order  Hankel  function  of  the  first  kind, 

(5.1)  Hn(z)  =  Jn(z )  +  iYn(z). 


5.1.  Functions  which  are  both  oscillatory  and  singular.  Classical  quadrature 
techniques,  like  Gaussian  quadratures,  perform  poorly  when  the  functions  to  be  inte¬ 
grated  exhibit  more  than  one  kind  of  oscillatory  or  singular  behavior.  In  this  example, 
we  present  quadrature  rules  for  functions  [0, 1]  — »  M  of  the  form 

n  /  m  \ 

<5-2>  E  E  dij  cos  (Pjx)  +  bij  sin  (fax)  I  xai, 

i= 1  \j= 1  / 
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where  the  cq  are  arbitrary  real  numbers  on  the  interval  [—.6,1.0],  and  the  are 
arbitrary  real  numbers  on  an  interval  [0,6]. 

To  construct  such  quadratures,  we  first  chose  fairly  large  integers  m  and  n,  and 
then  construct  n  Legendre  nodes  cq  on  the  interval  [—.6, 1.0]  and  m  Legendre  nodes 
/3j  on  the  interval  [0,6],  We  then  take  all  functions  of  the  form 

(5.3)  xai  cos  (fox)  and  xai  sin  (/3jx) 


Xi 

Wi 

Xi 

Wi 

0.1229785612931944E-21 

0.3552652759524009E-16 

0.1686329315315300E-12 

0.8100212080727273E-10 

0.8480396519577635E-08 

0.2999419813536928E-06 

0.4775750160241328E-05 

0.4182340572497058E-04 

0.2326791577716382E-03 

0.9141444687391752E-03 

0.2744584412742149E-02 

0.6682130398844386E-02 

0. 1379497029777539E-01 
0.2498028781481084E-01 
0.4073593719986401E-01 
0.6109753716986285E-01 
0.8573787689934898E-01 
0.1141325013227837E+00 

0. 1457010079575668E+00 
0.1798914175653836E+00 
0.2162162828182287E+00 
0.2542599995624801E+00 
0.2936726106825575E+00 
0.3341588668019299E+00 
0.3754666793859390E+00 
0.4173765534810784E+00 
0.4596923705410667E+00 
0.5022333728490615E+00 
0.5448270050049624E+00 
0.5873022022823650E+00 
0.6294826883836380E+00 
0.6711798141862015E+00 
0.7121844171914087E+00 
0.7522571077315104E+00 
0.7911163140906269E+00 
0.8284234137301751E+00 
0.8637645183766196E+00 
0.8966293685769854E+00 
0.9263901328780540E+00 
0.9522879806764252E+00 
0.9734441058145912E+00 
0.9889218194127054E+00 
0.9978661596654697E+00 

0.1975656354324724E- 
0.3567391763467614E- 
0.1205245279219683E- 
0.4310646162233048E- 
0.3435473776420015E- 
0.9380681263322151E- 
0.1165688301296981E- 
0.8041357440032319E- 
0.3553400421394265E- 
0.1117568231732080E- 
0.2706368315522836E- 
0.5351823901681276E- 
0.9027505156251501E- 
0. 1342713453729264E- 
0. 18088162261 18706E- 
0.2257750666507361E- 
0.2661334339415909E- 
0.3007780005914209E- 
0.3296666645566729E- 
0.3533282633911356E- 
0.3724775994063266E- 
0.3878153866772767E- 
0.3999467106048757E- 
0.4093600316907737E- 
0.4164319265743018E- 
0.4214400183613144E- 

0.4245762431761676E- 
0.4259573188716399E- 
0.4256313223505469E- 
0.4235800106227075E- 
0.4197166388463288E- 
0.4138788759007980E- 
0.4058161770611706E- 
0.3951708131079590E- 
0.3814519845080301E- 
0.3640037679024763E- 
0.3419715966891961E- 
0.3142816445555609E- 
0.2796677541655260E- 
0.2368149162837343E- 
0. 1847233461780559E- 
0.1233680051765303E- 
0.5456702457360249E- 

200 

L50 

L10 

)90 

)70 

)60 

)  10 

)  10 

)30 

)20 

)20 

)20 

)20 

)10 

)10 

)10 

)10 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)D. 

)29. 

7962322346848315E-22 
1830884530860187E-16 
7683378387780577E-13 
3548822655627213E-10 
3820064929774288E-08 
1453085581405940E-06 
2542069046177732E-05 
2452467963854058E-04 
1490989514657686E-03 
6325800496319742E-03 
2027015096036673E-02 
5213355986400703E-02 
1127123616180105E-01 
2121486403535344E-01 
3571945099051861E-01 
5498138040746230E-01 
7876554840715572E-01 
L065681257681138E+00 
L377839750139660E+00 
L718183483671372E+00 
2081399215101480E+00 
2462954187930989E+00 
2859051  752770733E+00 
S266515434452085E+00 
S682659628410250E+00 
U05169996086324E+00 
!531999641092431E+00 
1961279927664242E+00 
5391242133546008E+00 
5820145309561993E+00 
i246205472012950E+00 
566752101 1409383E+00 
TO81988727858907E+00 
(’487204161980902E+00 
r880339071872188E+00 
5257988675123893E+00 
5615983375780399E+00 
5949168202476468E+00 
5251176490208435E+00 
5514276420773934E+00 
S729462427273851E+00 
5887074670174289E+00 
S978240533043199E+00 

0.1255518558516502E-20 

0.1807251002083772E-15 

0.5430995592381254E-12 

0. 1887682739748173E-09 

0. 1568277392295615E-07 
0.4671196468523313E-06 
0.6453818186891716E-05 
0.4945575220250773E-04 
0.2400742808398220E-03 
0.8181326443575235E-03 
0.2119810413961939E-02 
0.4437537936998829E-02 
0.785221 1174980829E-02 
0.1215005865408980E-01 

0. 1689083767930460E-01 
0.2159081784351869E-01 
0.2589061574717723E-01 
0.2961 183182803991E-01 
0.3271991578613246E-01 
0.3526008081919886E-01 
0.3730762085807319E-01 
0.3894030217190292E-01 
0.4022648983646625E-01 
0.4122168685192172E-01 

0.4196875223655011E-01 

0.4249935402437477E-01 

0.4283556050595227E-01 

0.4299113717974382E-01 

0.4297240498182741E-01 

0.4277862015965356E-01 

0.4240185727864892E-01 

0.4182636166562592E-01 

0.4102730967187124E-01 

0.3996889290357553E-01 

0.3860165671827315E-01 

0.3685914346939630E-01 

0.3465427456539984E-01 

0. 318768749151 1273E-01 
0.2839583726541435E-01 
0.2407310551022889E-01 

0. 1880072022809006E-01 
0.1256986040620092E-01 
0.5563859224973946E-02 

TABLE  3.  Quadrature  formulas  for  the  functions  of  the  form  (5.2)  with 
a  G  [—.6, 1.0]  and  6  =  100.  The  43-point  rule  on  the  left  was  generated 
with  the  QR  variant  of  the  algorithm,  while  the  43-point  rule  on  the 
right  was  generated  with  the  SVD  variant.  Both  achieve  full  double 
precision  accuracy. 
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as  the  input  functions  <f>i  to  the  algorithm  of  Section  5.  It  is  clear  the  for  sufficiently 
large  n  and  m,  the  obtained  quadrature  will  work  for  all  functions  of  the  form  (5.2). 

In  Tables  1,  2,  and  3,  we  list  the  quadrature  nodes  and  weights  for  b  =  20,  50,  and 
100.  In  each  case,  two  quadrature  rules  were  computed,  one  using  the  QR  decompo¬ 
sition  variant  of  the  algorithm  and  one  using  the  SVD  variant.  The  parameters  m 
and  n  were  chosen  to  be  m  =  900  and  n  =  100.  In  Table  4,  we  report  the  time  spent 
computing  each  of  these  quadratures  as  well  as  the  number  of  quadrature  nodes  re¬ 
quired.  These  computations  were  performed  in  extended  precision  arithmetic  in  order 
to  ensure  double  precision  accuracy  for  the  resulting  quadrature  rule. 

We  also  computed  single  precision  (10-8)  quadrature  rules  for  b  =  20,  50,  and  100. 
These  computations  were  performed  in  double  precision  arithmetic  and  the  results 
are  reported  in  Table  5. 


b  = 

QR 

20 

SVD 

b  = 

QR 

50 

SVD 

b  = 

QR 

100 

SVD 

tcheb 

208.11 

213.21 

447.01 

457.93 

894.71 

933.67 

tnewt 

196.23 

195.03 

509.25 

634.77 

1277.34 

1315.04 

ttot 

404.94 

408.00 

956.14 

1091.22 

2171.12 

2248.01 

#  nodes 

26 

26 

33 

33 

43 

43 

Table  4.  Execution  times  (in  seconds)  for  the  computation  of  the 
double  precision  accuracy  quadratures  of  Example  5.1,  as  well  as  the 
number  of  quadrature  nodes.  Results  for  both  the  SVD  and  QR  vari¬ 
ants  of  the  algorithm  are  provided.  The  CPU  time  taken  by  the  first 
two  stages  of  the  algorithm  is  reported  as  tcheb  while  the  CPU  time  for 
Stage  3  is  reported  as  tnewt. 


b  = 

QR 

20 

SVD 

b  = 

QR 

50 

SVD 

b  = 

QR 

100 

SVD 

tcheb 

1.15 

1.18 

1.88 

1.94 

3.28 

3.26 

tnewt 

6.78 

6.43 

9.33 

9.86 

15.0 

22.3 

ttot 

7.93 

7.61 

11.2 

11.8 

18.2 

25.5 

#  nodes 

15 

15 

21 

21 

30 

30 

Table  5.  Execution  times  (in  seconds)  for  the  computation  of  the 
single  precision  accuracy  quadratures  of  Example  5.1,  as  well  as  the 
number  of  quadrature  nodes.  Results  for  both  the  SVD  and  QR  vari¬ 
ants  of  the  algorithm  are  provided.  The  CPU  time  taken  by  the  first 
two  stages  of  the  algorithm  is  reported  as  tcheb  while  the  CPU  time  for 
Stage  3  is  reported  as  tnewt. 
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5.2.  Plane  wave  expansions.  The  Green’s  function  for  the  Helmholtz  equation  in 
M3  satisfies  the  following  identity,  valid  for  z  >  0: 

iujr  poo  - - _  p 

(5.4)  —  =  /  e-2V«2-2J0(ev/^+^)^^=^ 

r  Jo  sje  -  J1 

where  r  =  \J  x2  +  y2  +  z2.  This  formula  can  be  derived  by  applying  the  Fourier 
Inversion  Theorem  followed  by  contour  integration.  In  [9],  a  scheme  for  accelerating 
fast  multipole  methods  for  the  Helmholtz  equations  at  low  frequency  was  introduced. 
It  operates  via  discretizations  of  formula  (5.4)  which  hold  for  x,  y,  and  z  satisfying 


L  <  z  <  4L 

(5.5)  — 4L  <  x,y  <  AL. 

The  appropriate  quadrature  rule  ostensively  depends  on  the  two  parameters  to  and 
L.  In  fact,  by  making  the  substitutions 

A  =  y  and  u  =  Luo 

Ju 

in  (5.4),  we  obtain  the  equivalent  representation 


(5.6) 


=  L 


e~{zL)'J^M\V(xL)2  +  ( yL )2) 


A 


:d\, 


'  i On 


which  shows  that  up  to  rescaling  factors,  the  quadrature  rule  depends  only  on  the 
product  uL ,  which  is  the  size  of  the  box  (5.5)  in  wavelengths.  Moreover,  since  Jq{z) 
satisfies  the  well-known  identity 


(5.7) 


J0(z)  =  —  cos(zsm(9))d9, 


7T 


which  can  be  found  as  formula  (9.1.18)  in  [1],  in  order  to  generate  a  discretization  of 
(5.4)  which  holds  for 

a  <  uL  <  6, 

it  suffices  to  take  as  input  to  the  algorithm  of  Section  5  functions  of  the  form 

„-W5J=5iA<»s(l/A) 


(5.8) 

where  x,  y ,  and  to  are  allowed  to  vary  as 


v/A^ 


( jj 


2  ’ 


1  <  x  <  4 

(5.9)  0  <  y  <  4 bV2, 

a  <  u  <  b. 


The  algorithm  of  this  paper  was  applied  to  functions  of  the  form  (5.8)  in  order  to 
generate  discretizations  of  formula  (5.4)  of  the  form 

picjr  _ _  c 

(5,10)  Ji - Wj 

r  i 

for  boxes  of  varying  sizes.  The  Stage  1  and  Stage  2  discretization  steps  were  per¬ 
formed  twice  for  each  choice  of  ujL  and  variant  (QR  or  SVD)  of  the  algorithm,  once 
using  extended  precision  arithmetic  and  once  using  double  precision  arithmetic.  The 
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c oL 

Expansion  order 

e  = 

QR 

10"3 

SVD 

e  = 

QR 

10“6 

SVD 

e  =  10"9 

QR  SVD 

e  = 

QR 

io-15 

SVD 

.5 

13 

13 

30 

29 

44 

44 

72 

72 

5 

18 

19 

33 

33 

46 

46 

80 

79 

10 

22 

24 

36 

35 

48 

49 

84 

84 

25 

44 

44 

53 

52 

65 

63 

94 

94 

Table  6.  Orders  of  the  plane  wave  expansions  of  Example  5.2  as  a 
function  of  box  size  and  accuracy. 


Stage  3  algorithm  was  then  used  repeatedly  to  generate  quadratures  of  varying  ac¬ 
curacies  for  each  of  the  boxes.  The  number  of  terms  in  the  plane  wave  expansion 
is  given  in  Table  6  as  a  function  of  the  size  of  the  box  in  wavelengths,  the  required 
accuracy,  and  which  variant  (QR  or  SVD)  of  the  algorithm  is  used.  Here  accuracy  is 
measured  as  the  largest  absolute  error  occurring  in  formula  (5.10). 

Table  7  reports  the  total  time  taken  by  the  Stage  1  and  Stage  2  computations. 
Finally,  Table  8  gives  the  time  taken  by  the  Stage  3  procedure.  Note  that  all  Stage 
3  computations  were  performed  in  double  precision  arithmetic,  except  for  the  quad¬ 
rature  rules  with  accuracy  10-15;  those  computations  were  performed  in  extended 
precision. 


lvL 

=  .5 

ojL 

=  5 

uL 

=  10 

uL 

=  25 

QR 

SVD 

QR 

SVD 

QR 

SVD 

QR 

SVD 

t  double 

8.84 

9.03 

8.23 

8.60 

7.17 

8.10 

9.14 

9.84 

tquad 

1858.96 

1954.31 

2046.89 

2192.72 

2079.71 

2227.13 

2486.69 

2644.03 

Table  7.  Time  (in  seconds)  taken  by  the  Stage  1  and  Stage  2  proce¬ 
dures  for  quadrature  rules  of  Example  5.2;  tdoubie  gives  the  CPU  time 
for  the  double  precision  computations  and  tquad  gives  the  total  CPU 
time  for  the  extended  precision  computations. 


5.3.  Integral  representations  for  //-expansions.  Let  (p,  9)  denote  the  polar  co¬ 
ordinate  system  defined  by 

x  =  pcos(6)  and  y  =  psin(d). 

As  is  well  known  (see,  for  instance,  [18]),  if  a  function  (j)  :  R2  — >  C  satisfies  the 
Helmholtz  equation 

V20  +  ujcj)  =  0 

outside  of  a  disc  D  of  radius  R  centered  at  0  and  the  radiation  condition 

lim  (f>(tx)e~tkt^Vi  =  c 

t— XX 


31 


c oL  CPU  Time  (seconds) 


e  = 

QR 

10“3 

SVD 

e  = 

QR 

10“6 

SVD 

e  = 

QR 

10"9 

SVD 

e  = 

QR 

io-15 

SVD 

.5 

0.58 

0.61 

3.19 

2.81 

15.01 

10.33 

12055.45 

12398.62 

5 

0.80 

0.86 

3.95 

3.58 

12.01 

11.32 

18468.92 

17513.46 

10 

1.25 

1.31 

5.03 

4.59 

14.31 

14.21 

22018.44 

22662.14 

25 

7.61 

8.01 

17.22 

1.57 

38.32 

35.21 

34646.45 

36221.12 

Table  8.  CPU  time  (in  seconds)  required  to  perform  the  Stage  3 
computations  for  the  quadratures  of  Example  5.2. 


at  oo,  then  <fi{x)  can  be  uniquely  represented  outside  the  disc  D  via  the  77-expansion 

OO 

(5.11)  <f>(x)  =  Y  PnHn(up)eine. 


Moreover,  once  N  >  \u\R,  the  error  in  the  approximation 

N 

(5.12)  <f>{x)  ~  Y  PnHn(up)eine. 

n=—N 


for  |  X  I  =  Ri  >  R  decays  as  (Ri/R)n  (see,  for  instance,  [18]). 
The  formula 


(5.13) 


Hn{cup)eind 


(-l)n  f°°  f  i£  -  -  £2 

77  J-  oo  -  £2  V  u 


e^dt, 


valid  for  y  >  0,  can  be  obtained  via  manipulation  of  the  well-known  representation 

(5.14)  Hn(z)  =  [  eizcos{w)+inwdw, 

71  Jc 

where  C  is  a  properly  chosen  contour  in  the  complex  plane  (see  [5]  for  a  simple 
derivation  of  (5.14)).  In  this  example,  we  construct  quadratures  for  formula  (5.13) 
which  hold  for  x  and  y  satisfying 


2L  <  y  <  AL 

(5.15)  —  L  <  x  <  L, 

and  for  n  —  0, 1, ... ,  M.  Since 

H-n(z )  =  (-1  )nHn(z), 

which  is  formula  (9.1.6)  in  [1],  an  L-point  quadrature  rule  for  functions  of  this  form 
allows  us  to  approximate  an  H-expansion 

M 

il>  =  Y  anHn{up)ein6 
n=—M 


(5.16) 
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,.17) 

-0  ~ 

b j(u 

3= 1 

uL 

M=2 

M=4 

00 

e  =  10- 

“7  e  =  10-15 

O 

1 

7  e  =  10~15 

e  =  10- 

-7  e  =  10-15 

1 

23 

57 

27 

63 

32 

66 

5 

29 

60 

29 

64 

33 

68 

10 

33 

66 

35 

67 

36 

72 

15 

39 

70 

39 

73 

41 

79 

Table  9.  Orders  of  the  expansions  (5.17). 


If  we  make  the  substitutions 


(5.18) 


A  =  £/L  and  u>  =  Luj0 


into  formula  (5.14),  then  we  obtain  the  equivalent  representation 


(5.19) 


Hn(up)eind 


(_l)n  roo  ei\(xL)+HyL)y/J*-\*  /^  _  ^2  _  y2 
7T  i-oo  y/^O  -  A2  V  U° 


d\, 


which  shows  that  as  with  the  last  example,  the  appropriate  discretization  for  (5.13) 
depends  only  on  uL. 

The  QR  variant  of  the  algorithm  was  used  to  compute  quadrature  rules  for  function 
of  the  form  (5.16)  for  various  values  of  uL  and  M.  Table  9  gives  the  number  of  terms 
in  the  expansion  (5.17)  necessary  to  obtain  precisions  10-'  and  1CT15  for  boxes  of 
different  sizes  and  for  different  values  of  M.  Table  10  gives  the  time  required  to 
perform  the  Stage  1  and  Stage  2  computations  for  these  quadrature  rules,  while 
Table  11  gives  the  time  taken  by  the  Stage  3  computations.  Note  that,  as  usual, 
the  computations  for  quadrature  rules  with  accuracy  10~'  were  performed  in  double 
precision  arithmetic  while  those  for  quadrature  rules  with  10~15  were  performed  in 
extended  precision  arithmetic. 


6.  Generalizations  and  Conclusions 

We  have  presented  a  simple  and  robust  scheme  for  the  computation  of  efficient 
quadrature  rules  for  a  wide  class  of  functions.  We  demonstrated  this  algorithm  by 
generating  efficient  quadrature  rules  for  several  different  classes  of  functions,  including 
systems  of  functions  exhibiting  different  kinds  of  singular  and  oscillatory  behavior. 

We  close  with  a  number  of  conclusions  and  possible  generalizations  of  this  work: 

1.  The  results  of  this  paper  are  purely  experimental.  While  there  is  a  framework  for 
proving  (under  certain  conditions)  the  existence  of  generalized  Gaussian  quadra¬ 
tures  (see  [11,  12,  15,  16,  13]),  it  does  not  apply  to  many  of  the  examples  of  this 
paper.  Indeed,  the  numerical  experiments  of  this  paper  and  those  of  [14,  3,  20] 
seem  to  indicate  that  such  quadratures  exist  under  very  general  conditions. 
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uL 

M 

=2 

M=4 

M= 

=8 

e  =  10“7 

e  =  10~15  e  =  10~7 

e  =  10~15 

e  =  10“7 

e  =  10”15 

1 

18.50 

1999.02 

24.16 

2796.14 

45.12 

4274.44 

2 

24.65 

2872.72 

32.96 

3811.83 

55.87 

5863.14 

5 

29.47 

3724.73 

41.62 

4693.70 

66.33 

7762.44 

15 

37.92 

5234.68 

50.98 

6359.60 

78.56 

8702.26 

Table  10.  Time  (in  seconds)  required  by  the  Stage  1  and  Stage  2 
computations  for  the  quadrature  rules  of  Example  5.3. 

ujL 

M=2 

M 

=4 

00 

(T\ 

o 

1 

e  =  10~15 

Cb 

O 

1 

-i 

e  =  10"15 

Cb 

o 

1 

-4 

e  =  10~15 

1 

1.01 

4238.30 

1.55 

6302.25 

2.85 

8411.11 

5 

2.08 

5694.78 

2.17 

6011.20 

3.26 

8302.36 

10 

3.38 

7604.56 

3.76 

6540.80 

4.36 

9221.55 

15 

6.39 

9881.15 

6.86 

10861.93 

6.94 

12322.11 

Table  11.  Time  (in  seconds)  taken  by  the  Stage  3  computations  for 
the  quadrature  rules  of  Example  5.3. 


2.  Stage  3  of  the  algorithm  of  Section  4,  wherein  quadrature  nodes  are  eliminated  one- 
by-one,  is  applicable  to  a  wide  range  of  problems.  The  method  applies  wherever 
a  sparse  solution  for  an  underdetermined  nonlinear  system  of  equations  is  being 
sought. 

3.  The  quadrature  weights  generated  by  the  algorithm  of  this  paper  are  generally, 
but  not  necessarily,  positive.  In  many  applications,  positive  quadrature  weights 
are  desirable.  A  modification  of  the  algorithm  to  ensure  that  the  weights  of  the 
resulting  quadrature  formulae  are  positive  might  prove  useful. 

The  Chebyshev  quadrature  procedure  of  Stage  2  of  the  algorithm  could  be  mod¬ 
ified  to  produce  positive  weights  by  replacing  the  least  squares  minimization  prob¬ 
lem  with  a  linear  program.  That  program  can  be  solved  via  the  simplex  method, 
which  would  result  in  a  Chebyshev  quadrature.  The  Stage  3  Newton  iterations 
could  be  modified  in  a  number  of  ways  (e.g.  barrier  methods)  to  ensure,  or  at  least 
encourage,  positive  quadrature  weights. 

4.  Although  there  are  several  difficulties  that  must  be  overcome,  there  is  no  funda¬ 
mental  barrier  to  generalizing  the  procedure  of  this  paper  to  higher  dimensions. 
A  procedure  for  generating  efficient  quadrature  rules  for  collections  of  functions 
defined  on  domains  in  M2  would  have  many  applications  in  numerical  analysis. 
This  topic  is  being  vigorously  investigated  by  the  authors. 
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5.  The  scheme  of  this  paper  has  obvious  application  to  the  interpolation  of  functions. 
In  particular,  the  observations  of  Subsection  2.3  and  the  algorithm  of  this  pa¬ 
per  allow  for  the  construction  of  stable,  efficient  interpolation  formulas  for  highly 
singular  and  oscillatory  functions. 

6.  There  are  a  number  of  applications  of  this  work  to  the  discretization  of  integral 
equations.  The  ability  to  produce  efficient  quadratures  for  broad  classes  of  col¬ 
lections  of  functions  is  by  itself  useful  in  the  discretization  of  integral  equations, 
and  the  general  framework  of  this  paper  for  “downsampling”  a  quadrature  for¬ 
mula  should  allow  for  the  reduction  of  the  complexity  of  discretizations  of  integral 
equations. 
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